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ABSTRACT 


The  research  deals  with  the  problem  of  apsferity-^excited  thermomechanical  field  in  a 
medium  with  a  surface  layer  and  a  near  surface  void  defect.  The  thermomechanical  field 
governs  the  mode  of  cracking,  which  leads  to  failure  in  the  wear  surface.  The  presence 
and  location  of  the  void  defect  is  most  critical.  This  investigation  obtained  the  solutions  for 
the  temperature  distribution  and  the  stress  state  in  a  layered  medium  with  a  rectangular 
cavity.  This  temperarture  distribution  and  stress  state  result  when  the  solid  medium  is 
subjected  to  Coulomb  frictional  loading  from  an  asperity  moving  at  a  moderately  high 
speed  (of  approximately  10-15  m/s).  In  the  analysis,  the  coated  medium  was  represented 
by  a  solid  half  space,  with  a  thin  top  surface-layer  of  solid  wear  material.  The  cavity 
defect  required  a  mathematical  model  in  terms  of  the  material  coordinates.  The  corre¬ 
sponding  governing  differential  equations  were  time-explicit  and  transient.  A  general  finite 
difference  formulation  was  developed  to  calculate  both  the  temperature  and  the  stress 
fields.  The  energy  balance  method  was  applied  at  the  corners  of  the  rectangular  cavity  to 
resolve  the  problem  of  singularities  in  the  temperature  field.  JThe  stress  singularity  at  each 
corner  was  represented  by  a  special  element  that  was  introduced  representing  the  behavior 
of  the  known  stress  singularity  at  the  corner  and  its  vicinity.  The  general  equation  of  the 
stress  field,  including  the  dynamic  term,  is  of  the  regular  perturbation  type.  The  small 
order  dynamic  term  is  demonstrated  to  be  a  higher  order  effect  by  perturbation  method, 
thus  negligible.  Numerical  solutions  were  carried  out  for  the  zeroth  order  approximation 
and  the  case  of  uniform  asperity  pressure  distribution. 

It  was  shown  that,  at  moderately  high  asperity  speed,  the  thermal  stress  effect  domi¬ 
nates  the  combined  thermo-mechanical  stress  field,  which  eventually  leads  to  failure  in  the 
no-cavity  case.  When  a  defect,  such  as  a  cavity,  exists,  the  stress  state  that  determines  the 
failure  phenomenon  is  much  more  severe  and  can  be  quantified  depending  on  the  location 
of  the  cavity.  These  results  are  determined  through  a  numerical  computation  based  on  the 
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material  properties  of  Stellite  111.  However,  the  parametric  effects  of  material  variations  in 
the  coating  and  the  substrate,  including  changes  in  both  thermal  and  mechanical  properties, 
were  also  considered.  The  study  of  the  cavity  location  also  established  the  existence  of  a 
critical  cavity  location  for  cracking  by  cohesive  failure.  This  location  is  defined  by  the 
critical  ligament  thickness  (thickness  between  the  wear  surface  and  the  top  edge  of  the 
cavity),  at  which  the  cavity-influenced  thermal  tensile  stress  reaches  a  maximum.  This 
thickness  is  important  to  designers  when  cavities  at  coating/substrate  interfaces  are  either 
unavoidable  or  too  expensive  to  control  in  fabrication. 
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INTRODUCTION 


1.1  Statement  of  Problem 

This  Investigation  studies  the  thermomechanical  cracking  in  a 
coated  Medium  with  a  near  surface  cavity.  Such  cavity  generally  occurs 
in  the  neighborhood  of  the  coating/substrate  interface,  as  a  result  of 
either  Inclusion  or  poor  bonding  during  the  coating  process.  A  typical 
geometry  of  the  cavity  can  be  shown  as  in  Figure  1.1.  To  facilitate 
the  analysis,  this  research  will  first  study  the  effect  of  a 
rectangular  cavity.  The  general  failure  mechanism  is  caused  by  the 
frictional  excitation  of  a  moderately  high  speed  asperity  traversing 
over  a  coated  surface.  The  understanding  of  this  failure  process  will 
improve  the  design  of  the  modified  wear  surface  by  alleviating  the 
problem  of  friction  cracking  or  delamination. 

When  two  fiat  solids,  which  are  placed  in  contact  under  heavy 
loads,  slide  relative  to  each  other,  the  nominal  design  pressure 
between  the  mating  surfaces  is  based  upon  the  nominal  design  contact 
area.  When  the  contact  pressure  is  evenly  distributed  according  to 
design,  the  service  life  of  the  solids  is  not  a  serious  problem,  even 
at  a  high  rubbing  speed.  However,  at  high  operating  speed,  the  real 
contact  area  can  be  reduced  by  several  orders  of  magnitude.  As  a 
result,  a  low  design  pressure  may  result  in  a  very  high  Interfacial 
pressure,  thus  a  very  high  dry  frictional  force  in  the  real  contact 
area.  Kennedy  111  shoved  that  the  size  of  the  real  contact  area 
depends  on  operating  speed  and  material  parameters  such  as  thermal 


conductivity  and  wear  resistance.  He  concluded  that  decreasing  thermal 
conductivity,  increasing  wear  resistance  and  Increasing  operating  speed 
will  reduce  the  real  contact  area.  It  was  shown  that,  for  a 
conservative  areal  ratio  (contact  area/nominal  area>  of  10“3  (Burton 
121  considered  10”  4  as  a  possible  areal  ratio),  a  low  design  pressure 
of  240  kPa  (35  psi)  could  result  in  a  240  MPa  (35,000  psi)  local 
pressure  in  the  contact  zone.  The  high  friction  would  generate  locally 
an  extremely  high  temperature,  which  was  called  "flash  temperature"  by 
Archard  131.  The  local  contact  area  is  called  "red  banding"  or  "hot 
spot"  C41,  which  has  been  experimentally  demonstrated.  In  severe  cases 
the  temperature  can  be  extremely  high,  leading  to  cracking  of  the 
surface  151.  This  phenomenon  is  called  "heat  checking"  or 
"thermocracking"  163.  It  is  frequently  seen  in  seal  rings,  brakes,  and 
rail-wheels  [7,8,9,10,11,121  as  shown  in  Figures  1.2  and  1.3.  In 
general,  it  was  observed  that  numerous  radial  cracks  developed 
perpendicular  to  the  sliding  direction  and  almost  periodically  along 
the  circumference.  In  order  to  understand  these  failures,  in  recent 
years,  there  has  been  increased  emphasis  in  finding  a  solution  of 
failure  control,  both  experimentally  and  analytically. 

1.2  Related  Investigation  in  Progress 

The  phenomenon  of  high  temperature  "hot  spot"  was  observed  in  the 
experiments  by  Archard  (33.  A  general  survey  of  the  problem  of 
cracking  through  the  development  of  a  frictional  hot  spot  was  discussed 
by  Burton  143.  Proof  of  the  existence  of  hot  patches  of  solid-to-solid 
contact  was  obtained  experimentally  by  Banner Jee  and  Burton  1133  in  the 
case  of  metallic  rings  rotating  against  a  non-metal lie  disk  and,  more 
recently,  in  actual  operating  face  seals  by  Kennedy  C141.  The  latter 


Figure  1.3  Thermal  cracks  on  the  brake  shoe 
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study  Bade  use  of  a  new  contact  probe  which  enables  the  monitoring  of 
contact  patch  sizes  and  locations  in  ring-on-ring  or  ring-on-disk 
configurations.  It  was  proven  to  be  quite  effective  in  determining  the 
geometry  and  movement  of  contact  patches  in  dry  operation  of  mechanical 
face  seals.  In  his  earlier  experiment  L153,  Kennedy  used  a  carbon  ring 
against  a  metallic  mating  ring  made  from  440  C  stainless  steel, 
beryllium  copper  or  52100  bearing  steel  under  both  dry  and  liquid 
lubricated  conditions.  In  these  experiments,  the  existence  of  distinct 
spot  asperities  on  the  metallic  ring  was  also  observed.  It  was  found 
that  the  spots  tend  to  remain  stationary  with  respect  to  the  metallic 
mating  ring  of  the  seal,  whether  that  ring  is  stationary  or  rotating. 
However,  other  investigations  have  shown  hot  patches  moving  relative  to 
the  mating  ring  and  stationary  on  the  primary  ring  (13,161.  Burton 
[17]  also  reported  that,  for  an  aluminum  ring  sliding  on  a  glass  disc, 
the  hot  spot  precessed  at  a  much  lower  speed  than  the  rubbing  speed. 

The  uncertainty  of  this  observed  discrepancy  on  the  speed  of  the  moving 
asperities  remains,  but  there  Is  no  doubt  about  the  existence  of  the 
moving  asperities  due  to  thermoelastlc  Instability  on  mechanical  face 
seals.  Several  analytical  studies  of  the  failure  due  to  the  existence 
of  the  moving  asperities  have  been  developed.  Surface  displacements, 
temperature  field  and  stress  state  of  a  convective  elastic  half  space 
under  an  arbitrarily  distributed  fast-moving  line  heat  source  were 
obtained,  using  Integral  transform  techniques,  by  Ling  et  al  [18,19,201 
and  How  and  Cheng  C 21 3 .  Kilaparti  and  Burton  C221  have  developed  an 
exact  Fourier  series  solution  for  a  periodic  strip  heat  input.  Their 
series  is  rather  unwieldly,  but,  at  large  Peclet  number  (R=Va/x),  it 
reduces  to  a  form  (231  that  is  simpler  than  that  of  Ling  and  Mow  (181. 


Recently,  Barber  £241  employed  the  Green's  function  for  the  problem  of 
Kllaparti  and  Burton,  and  obtained  the  thermoelastlc  displacements  and 
stresses  due  to  a  heat  source  moving  over  the  surface  of  a  half  plane. 

A  finite  element  analysis  was  developed  by  Kennedy  C 25 ]  to  study  the 
surface  temperatures  resulting  from  frictional  heating  in  sliding 
systems.  He  also  applied  finite  element  techniques  to  study  the 
stresses  in  the  mechanical  face  seals  161  and  showed  that  the  dominant 
stresses  in  the  seal  components  are  thermal  stresses.  The  surface 
stress  component  (parallel  to  the  surface)  resulting  from  a  periodic 
row  of  moving  hot  patches,  with  width  2a  each,  and  a  spacing  of  2m  was 
Investigated  by  Tseng  and  Burton  1261.  They  concluded  that  the  tensile 
stress  would  appear  instantaneously  with  each  passage  of  the  heat 
source.  Two-  dimensional  models  of  heat  checking  in  the  contact  zone 
of  a  face  seal  were  presented  by  Ju  and  Huang  [27].  Because  of  the 
three-dimensional  aspect  of  those  observed  "hot  spots",  Ju  and  Huang 
reformulated  the  problem  in  three-dimensional  theory  of 
thermoelasticity  128,29,301.  The  investigation  concluded  that  the 
highest  tensile  stress  occurs,  for  an  asperity  speed  of  10-15  m/s  <400- 
600  ln/s),  at  a  depth  of  the  order  of  one-tenth  the  asperity  size. 

This  depth  defines  the  critical  depth  of  the  material.  The  physical 
depth  is  therefore  50-100  pm.  At  such  a  asperity  speed,  the  stresses 
from  the  thermal  effect  of  the  asperity  friction  are  an  order  of 
magnitude  larger  than  those  from  its  mechanical  traction  effect.  Ju 
and  Huang  C 31 J  also  demonstrated  that,  when  asperities  excite  the 
surface  periodically  in  close  Intervals  (a  numerical  example  used  a 
spacing  of  twelve  asperity  size),  the  thermomechanical  effects 
accumulate,  yet  tending  to  a  limiting  magnitude,  even  though  the 


mechanical  stress  dissipates  with  no  residue  effect.  The  cumulative 
effect  definitely  depends  on  the  interval  of  periodic  excitations.  At 
a  relatively  large  interval  of  approximately  1000  asperity  size,  no 
cumulative  effect  is  evident. 

For  improvement  of  the  wear  property  of  the  surface,  recent  effort 
has  been  directed  toward  surface  modifications.  Research  to  understand 
the  behavior  of  coated  surfaces  under  asperity  excitation,  hence,  has 
gained  importance.  Ju  and  Chen  t 32, 333  first  solved  for  the  case  of  a 
moderately  thick  coating  (thickness  of  the  order  of  the  asperity  size). 
Later  Ju  and  Liu  C341  extended  the  general  formulation  of  [32,331  to 
study  the  thickness  effect  of  the  coating  layer  for  various  mechanical 
and  thermal  Impedance  matchings  between  the  surface  coating  layer  and 
the  substrate.  It  is  concluded  by  Ju  et  al  that:  (i)  a  stiff  surface 
layer  would  result  in  higher  thermal  stress;  (11)  the  stress  state  in 
layered  media  is  Influenced  by  the  layer  thickness,  reaching  a  worst 
state  when  the  coating  layer  thickness  is  in  the  neighborhood  of  the 
critical  depth;  (ill)  a  substrate  of  lower  thermal  expansion 
coefficient,  higher  Young's  modulus,  higher  thermal  conductivity  and 
capacity  will  result  in  lower  stresses  in  the  coating  layer;  (iv)  for 
the  thin  coating  layer,  the  shearing  stress  at  the  coating/substrate 
Interface  is  by  no  means  trivial,  depending  again  on  the  surface 
coating  thickness.  The  interface  shear  reaches  a  maximum  when  the 
coating  thickness  is  in  the  neighborhood  of  the  thermal  layer.  These 
results  are  important  for  designing  the  bonding  of  the  surface  coating. 

In  the  previous  work  on  the  moving  asperity  problem,  the  analyses 
dealt  with  basically  uniform  solid  media;  that  Is,  the  material  and 
asperity  properties  are  Invariant  in  the  direction  of  the  asperity 


motion.  In  such  cases,  since  the  time  effect  can  be  rendered  implicit 
in  the  Fourier  and  the  Navier  equations  by  using  a  coordinate  system 
fixed  to  the  traversing  asperity  (called  the  convective  coordinate 
system),  the  resulting  solutions  are  steady-state.  However,  when  the 
material  has  a  cavity,  uniformity  in  the  direction  of  the  asperity 
motion  no  longer  exists.  Consequently,  a  coordinate  system  fixed 
either  to  the  cavity  or  to  the  material  (referred  to  as  the  material 
coordinates)  must  be  employed.  The  governing  equations  and  their 
solutions,  therefore,  are  transient.  The  present  investigation  not 
only  obtains  the  temperature  field  solutions  but  also  analyzes  the 
stress  field  caused  by  the  input  of  a  moving  heat  source.  In  this 
study,  since  the  Fourier  and  the  thermoelastic  Navier’ s  equations  in 
the  material  coordinates  are  time  explicit,  the  finite  difference 
method  is  considered  more  appropriate.  Although  a  specific  numerical 
solution  does  not  show  the  effects  of  parameters,  a  general  trend  of 
the  parameters  effects  can  be  obtained  with  adequate  numerical 
solutions  for  a  series  of  given  parametric  values. 

1.3  General  Theory 

The  phenomenon  of  thermomechanical  cracking,  as  observed  from 
experiments  and  operational  damages,  is  connected  with  relatively  hard 
materials;  such  as  cast  iron  and  Stellite  III.  Blau  1351  and  Ruff  and 
Blau  1363  demonstrated  experimentally  that  the  plastic  wear  and  surface 
shear  for  hard  wear  material  are  restricted  to  a  very  thin  surface 
layer  (about  4-8p).  Ju  et  al  [27,28,29,30,33,34]  also  proved  that  the 
critical  depth  is  at  a  depth  of  an  order  of  magnitude  larger  than 
plastic  depth.  Therefore,  the  linear  thermoelastic  theory  holds.  The 
basic  mathematical  formulation  of  uncoupled  thermoclastlcity  consists 
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of  the  following  equations: 


pV2u  ♦  <X+p)  grad  div  u  -  <3X+2p>  a  grad  T  =  pu  (1.1) 

and  kV2T  =  pci'  (1.2) 

where  T  and  u  are  temperature  and  displacement  fields,  respectivity , 
k  is  the  thermal  conductivity,  p  is  the  mass  density,  c  is  the 
specific  heat,  X,p  are  the  Lame  constants,  and  a  is  the  coefficient 
of  thermal  expansion.  The  coupling  term  is  negligible  except  for 
conditions  in  which  the  temperature  distributions  have  sharp 
variations  in  their  time  histories,  which  often  occurs  during  the 
propagation  of  thermoelastic  waves  in  the  aftermath  of  thermal  shocks 
£37,38,39,40,41,42].  For  the  current  problem,  since  the  asperity  speed 
under  consideration  is  much  slower  than  the  elastic  wave  speed,  the 
uncoupled  thermoelastic  theory  is  applied. 

The  dynamic  effect  may  result  from  either  a  dynamic  loading  state 
or  a  non-steady  thermal  state  in  which  the  time  rate  of  temperature 
change  could  keep  up  with  the  stress  waves  in  the  material.  Duhamel 
£43]  stated  that  the  inertia  term  can  be  disregarded  if  the  time  rate 
of  change  of  temperature  is  slow  enough.  Parkus  £44]  showed  that  the 
significant  effect  from  the  inertia  term  can  arise  only  when  there  is 
an  Instantaneous  change  in  the  surface  temperature  or  in  the 
temperature  of  the  surrounding  medium.  In  fact,  the  dynamic  effect  is 
greatly  reduced  if  the  temperature  change  occurs  in  a  very  short,  but 
finite,  interval  of  time.  This  was  confirmed  by  Danilovskaya  £45,461, 
who  studied  the  dynamic  effect  due  to  a  thermal  shock  on  the  surface  of 


a  half-space  and  demonstrated  that  the  maximum  dynamic  stress  Is 
reduced  to  86Z  even  for  the  extremely  short  duration  heating  of  10”12 
seconds.  In  general,  under  usual  conditions  of  heat  exchange,  the  rate 
of  temperature  change  is  small  in  comparison  with  the  speed  of  sound  in 
the  material.  Thus,  at  any  instant,  the  thermal  stress  state  can  be 
determined  by  the  instantaneous  values  of  the  temperature  field. 

For  the  cavity  problem,  the  effect  of  the  dynamic  term  in  Equation 
(1.1)  will  be  studied  quantitatively  with  a  perturbation  method.  That 
is,  the  solution  to  Equation  (1.1)  can  be  expressed  in  an  asymptotic 
series.  Substituting  this  series  into  Equation  (1.1)  leads  to  a  set  of 
linear  equations  for  u.  Each  set  of  linear  equations  represents  a 
different  order  of  solution  of  the  asymptotic  series.  The  details  of 
the  perturbation  procedure  will  be  addressed  in  Chapter  4. 


ANALYTICAL  MODEL  AND  BASIC  EQUATIONS 


The  experiments  performed  by  Kennedy  147]  have  shown  that  contact 
between  two  flat  conforming  rings  is  concentrated  in  several  (1  to  5> 
patches,  with  a  few  small  solid-solid  contact  spots  occurring  within 
each  patch.  Each  contact  spot  is  identical  and  the  contacts  are 
equally  spaced  around  the  ring  circumference.  A  ring  could  therefore 
be  divided  Into  as  many  sections  as  the  number  of  contact  spots  and 
only  one  such  section  would  have  to  be  analyzed.  Kennedy  [15]  also 
proved  that  the  width  of  the  contact  spot  (asperity)  Is  about  0.1  to  1 
mm  <0.004  to  0.04  in.);  however,  the  size  of  a  typical  mating  ring  is 
several  orders  of  magnitude  larger  than  the  asperity  size.  Because  of 
this  size  difference  between  the  contact  area  and  the  mating  rings,  the 
analytical  model  is  represented  by  a  semi-infinite  body  with  a  thin 
coating  layer  and  a  rectangular  cavity  in  the  neighborhood  of  the 
coating  layer/substrate  Interface.  The  half  space  surface  is  subjected 
to  the  frictional  heating  of  a  moving  asperity  over  the  wear  surface 
(Figure  2.1),  and  the  material  coordinate  system  (fixed  to  the  cavity) 
is  used.  As  presented  in  Chapter  1,  the  linear  thermoelastic  theory 
applies  for  the  current  problem.  The  advantage  of  the  linear  theory  is 
the  application  of  the  superposition  principle,  which  allows  a 
separation  of  the  stress  field  to  a  contribution  of  the  mechanical  load 
of  the  pressure  and  friction  from  the  moving  asperity  and  another 
contribution  of  the  heat  input  from  the  rate  of  the  frictional  energy 
dissipation.  The  combined  effects  will  then  determine  the  possibility 
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of  fracture  Initiation.  The  governing  differential  equations  for  the 
temperature  and  the  stress  fields  are  the  Fourier  equation  and  the 
thermoelastic  Navier's  equation,  respectively. 

2.1  Temperature  Field 

The  governing  equation  for  the  temperature  field  is 


(2.1) 


where  T  is  the  temperature;  k  is  the  thermal  conductivity;  p  is  the 
mass  density;  c  is  the  specific  heat;  and  3  denotes  the  layered 
region.*  I  for  the  coating  surface,  II  for  the  substrate. 

The  temperature  field  T  must  satisfy  the  initial  condition. 


D 

TtXj.Xg.O)  =  0, 


(2.2) 


and  the  boundary  conditions: 

(i>  The  regularity  condition  holds  at  infinity  (Xj2+x22-*»>  , 

TP  =  0.  (2.3) 


(ii)  In  the  asperity  contact  surface  (c(t><  xt  <  c(t)+a,  x2=0),  the 
maximum  heat  input  through  the  boundary  is  the  rate  of  the  frictional 
energy 


-  k. 


3T 

3x2 


q  *  PfVp' (Xj ) , 


(2.4) 


where  is  the  Coulomb  coefficient  of  friction;  V  is  the  asperity 
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velocity;  p'(Xj)  is  the  pressure  over  the  contact  area;  q  is  the  heat 


flux  through  the  contact  area;  and  c(t)  is  the  distance  from  xt 


origin  to  the  leading  edge  of  the  asperity. 


<iii)  Outside  the  contact  surface,  xt  <  c(t>  or  xt  >  c<t)+a,  x2  =  0, 


the  convective  heat  loss,  being  of  small  order,  is  neglected  without 


loss  of  generality, 


3T 

5*1  =  °- 


(2.5) 


(iv)  At  the  coating  layer/substrate  Interface,  x2  =  H,  the  continuity 


conditions  hold 


T1  =  T11  , 


(2.6) 


3T1  at1* 


.  v  a  •  v  a 

kl  5*^  *  55T 


<  2. 7  > 


where  H  is  the  layer  thickness. 


<v>  Adiabatic  conditions  at  cavity  boundaries, 


5^  =  0,  at  -d  <  x,  <  d,  x2  =  L* 


(2.8) 


5x^~  =  °.  at  "d  i  X1  i  d’  X2  ■  L'+e 


(2.9) 


gj —  =  0,  at  Xj  =  d ,  L'  £,x^$.L'+e 


(2.10) 


IxT  =  0,  at  x.  *  -d,  L’  <  x2  <  L'+e 


(2.11) 


where  <1  is  the  half  width  of  the  cavity  and  e  is  the  depth  of  the  cavity. 
The  region  between  the  cavity  and  the  wear  surface,  which  is  important  in 
determining  the  magnitude  of  the  temperature  field  and  the  stress  state, 
shall  be  designated  descriptively  as  the  ligament  region.  The  distance 
between  the  surface  and  the  cavity  top  edge  is  therefore  the  ligament 
thickness  L'  (see  Figure  2.1). 

2.2  Mechanical  stress  field 

The  elastic  Navier's  equation  and  the  Hooke's  law  equation  are 
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i,J,k  =  1,2 
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(2.12) 
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a?  =  X  6  <  -is  )  +  v  (  -A  +  —1  ),  i,j,k  -1,2  (2.13) 
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where  ut  and  u2  are  the  displacements  in  xt  and  x2  direction, 
respectively;  d1J,d12,o22  are  stress  components;  X  and  p  are  the  Lame 
constants;  and  d  is  the  Kronecker  delta.  The  summation  convention  is 
used  for  all  repeated  indices  of  Roman  minuscules. 

The  mechanical  stress  field  is  initially  homogeneous. 

The  boundary  conditions  are: 

(1)  On  the  contact  surface  (c(t)  <  Xj  <  c(t)+a,  x2  =  0),  tractions 
are  prescribed  by 


pfp’ <xt ) , 


(2.14) 


(2.15) 


(ii)  Outside  the  contact  area,  the  region  given  by  xt  <  c(t)  or 
x1>c(t)+a,  x2  =  0,  the  surface  tractions  are  Identically  zero: 


(2.16) 


(2.17) 


(iii)  The  regularity  conditions  hold  at  infinity  (x12+x22-*«> , 


u*  =  0, 


(2.18) 


o*j=0. 


(2.19) 


(iv)  At  the  coating  layer/substrate  interface,  x2  =  H,  the 
continuity  conditions  are 


I  II 

u  =  u  , 
i  i 


I  II 
o=o 
12  12 


(2.20) 


(2.21) 
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(v)  The  cavity  boundary  is  traction  free;  that  is. 


=  °22  =  °*  at  -d  <  xt  <  d  ,  x2  =  L' 


(2.22) 


°i2  ~  °22  =  0»  at  -d  <  xt  <  d  ,  x2  =  L'+e 


(2.23) 


m 


°ii  =  °i2  =  0.  at  xt  =  d  ,  L'  <  x2  <  L’+e 


(2.24) 


Id 


(2.25) 


o^t  =  o^2  =  0,  at  Xj  =  -d  ,  L'  <  x2  <  L’+e 

2.3  Thermal  Stress  Field 

The  thermoelastic  Navicr's  equation  is 
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(2.28) 


and  the  Hooke's  law  equation  is 


3u; 


3uP  3uP 


«  .  (  :rJj-  )  ♦  u.  (  =-*-  ♦ 
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-  <3xp  +  2V  %  *tJ  <T-T0>  . 


(2.27) 


where  a  is  the  coefficient  of  thermal  expansion,  T  and  its  derivatives 
are  derived  from  the  temperature  field. 

The  initial  conditions  for  the  thermal  stress  field  are 


u*  ( x t , x2 , 0 )  =  0,  (2.28) 

°ij(xi ,x2,0)  =  0.  (2.29) 

The  boundary  conditions  are: 

(1)  The  surface,  xz  =  0,  is  traction  free,  i.e. 

o[2  =  0,  (2.30) 
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022  =  0. 


(2.31) 


<ii)  The  regularity  conditions  hold  at  infinity,  x^+Xj2-**,  i.c. 


uP  =  0, 


(2.32) 


(2.33) 


(ill)  Continuity  conditions  hold  at  the  interface,  x~  =  H,  i.e. 


t  II 
Ui  *  \  ’ 


(2.34) 


I  _  II 
a  =  o 
12  12 


(2.35) 


(iv>  The  cavity  boundary  is  traction  free,  i.e. 


0^2  =  °22  =  at  -d  <  <  d  ,  x2  =  L' 


(2.36) 


°i2  =  °22  =  0»  at  -d  <  xt  <  d  ,  x2  *  L’+e 


(2.37) 


=  o^2  =0,  at  Xj  -  d  ,  L'  <  x2  <  L*+e 


(2.38) 


o?t  =  o^2  =  0,  at  xt  =  -d  ,  L'  <  x2  <  L'*e 


(2.39) 


The  solution  techniques  and  the  numerical  results  of  the 


temperature  and  the  stress  fields  shall  be  given  in  Chapters  3  and  4 
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CHAPTER  3 

TEMPERATURE  SOLUTIONS 

Since  a  high  temperature  and  its  gradients  are  the  source  of  the 
high  thermal  stresses,  which  can  lead  to  the  thermocracking  of  the  wear 
medium.  Therefore,  it  is  of  primary  importance  that  a  temperature 
solution  is  available.  The  governing  equations  in  the  dimensionless 
form  are  as  follows: 

In  the  coating  layer,  0<n<D,  denoted  by  the  superscript  I, 


32*1  32*1 

ac2  + 


a#1 

an2  =  Ri 


(3.1) 


In  the  substrate  region,  D<n<®,  denoted  by  superscript  II, 


32*11  32*11 


as2 


an2 


a*11 

=  Rn  3t  * 


(3.2) 


where  ♦<Wk  r/q0a)  is  the  dimensionless  temperature;  ((,n)=(xt/a,  x2/a) 

are  the  dimensionless  coordinates  in  the  direction  opposed  to  the 

asperity  motion  and  the  depth  direction,  respectively  (as  shown  in 

Figure  2.1);  x=(Vt/a)  Is  the  dimensionless  time;  D=H/a  is  the 

dimensionless  coating  thickness;  R  =(Va/x_>  are  the  Peclet  numbers  in 

P  P 

the  coating  layer  (3=1)  and  in  the  substrate  ($=II>;  q0  is  the  average 
heat  flux  through  the  contact  area;  and  T,  K,  k,  a,  x|(  x2.  V,  t,  H  are 
the  same  as  defined  in  Chapter  2. 

3.1  Difference  Formulation 

Because  of  the  analytical  complexity  of  the  mathematical  model,  the 
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explicit  finite  difference  method  is  employed  to  solve  the  current 
problem  [48,49,50, SI , S2, S3, S4, SS, ] .  A  brief  discussion  of  the  finite 
difference  method  is  given  in  Appendix  I.  In  the  finite  difference 
method,  the  semi-infinite  body  is  replaced  by  a  sufficiently  large 
rectangular  region  (Figure  3,1).  A  central  difference  is  used  for  the 
space  derivatives,  and  a  two-point  forward  difference  is  used  for  the 
time  derivative  of  the  first  time  step,  then  a  three-point  forward 
difference  is  employed  for  the  following  time  steps.  The  reason  for 
using  the  three-point  forward  difference  after  the  first  time  step  is 
that  it  is  more  accurate  than  the  two-point  forward  difference.  But, 
for  the  first  time  step,  we  have  Information  only  on  one  previous  time 
line  (initial  condition),  and,  therefore,  only  the  two-point  difference 
formula  may  be  used. 

The  governing  differential  equations  in  the  difference  form  are: 

In  the  coating  layer 

♦I(i,J,n)  =  r1*I(i-l, J,n-1)  ♦  Cl  -  2(rt  ♦  r2)l*I(i, J,n-1)  + 

+  rl*I(i+l, J,n-1)  ♦  r2*I(i , J-l,n-l>  + 

♦  r2#J(i, j+l,n-l ) ,  n=l  (3.3a) 


and 


T  1  T  T 

*  =  r  [-♦<!, j,n-2>  ♦  2ra*A<i-lf J,n-1>  ♦ 


♦  4(1  -  r,  -  Tj)*^!, J,n-1>  +  2r j*1  ( i+l , j ,n-l )  ♦ 
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*  2r2*T(i, j-l,n-l>  +  2r2*r< i, j+l ,n-l ) 3  ,  n>l  <3.3b> 
In  the  substrate 

♦XI(i,j,n>  =  fl2r1*I1(i-l, J,n-1)  +  [1  -  ( 2Q2r,+2G2r2 > 1*IJ( i , j , n-l >+ 

♦  02r1*II(i+l, j,n-l)  ♦  02r 24x 1 < i , j-l , n-1 )  + 

♦  02r2*1I(i,3+l,n-i>  ,  n=l  (3.4a) 

and 

♦IT(i, J,n>  =  g  (-*lI(i, J,n~2>  ♦  20*rtfXI<i-l,jtn-l>  ♦ 

♦  4(1  -  Q2r  j  -  Q2r2>*11( i , j ,n-l >+2Q2r14II< i+1, j , n-1 )  ♦ 

+  202r2*II(i, j-l,n-l)+2n2r2#ir(i, j*l  n>2  (3.4b) 

where  r1*A'c/(Rj •  A^2 >  f  r2=A'c/(RIJ •  An2 > ,  °2=KI1/|CI*  *nd  denotes 

the  two  spatial  indices  and  the  time  step,  respectively.  For  the 
explicit  scheme,  the  time  step  Lx  must  satisfy  a  stability  criterion. 
The  most  commonly  used  method  of  stability  analysis  is  Von  Neumann's 
method  (48,50,54,551.  In  this  method,  a  finite  Fourier  series 
expansion  of  the  solution  to  a  model  equation  is  made,  and  the  decay  or 
amplification  of  each  mode  is  considered  separately  to  determine 
stability  or  instability,  as  we  now  demonstrate. 

Consider  first  the  difference  form  of  Equation  (3.1)  of  the 


coating  layer 


=  ♦I(i,J,n-l)  ♦  r,!*1  ( i-1, j ,n-l )  -  2*z(l.j,n-l>  ♦ 

+  ♦^i+l.J.n)]  +  r^^a.j-l.n-l)  -  2*r(i,J.n-l)  ♦ 

♦  ♦I(i, j+1  ,n-l)3  .  ( 3. 55) 


Each  Fourier  component  of  the  solution  is  written  as 


♦x<i, J,n> 


=  VneJK^<lAC>eJI:n<4Ar,> 


<3. 6a) 


where  Vn  is  the  amplitude  function  at  time-level  n  of  the  particular 


component  whose  wave  numbers  are  and  in  E  and  n  directions  and 
J=/-i  .  If  ©=K..AE  and  ©=K  An  we  obtain 

«.  n 


vneJ<iet>j^> 


(3.6b) 


Substituting  Equation  (3.6b>  into  Equation  (3.5)  gives 

vneJ(i6+j$)  =  vn-ieJ(ie+J©)  +  rilvn-ieJ[(i-i>e+J*l_2vn-ieJ(ie+^)  + 

♦  Vn_1eJ^ { 1  j  4  r2{Vn_lejCie+<  ^“i^]  _ 

-  2Vn‘1e<i®+^>  ♦  vn“teJlie+<^1>*3}  .  (3.7) 

Canceling  the  comsK>n  term  gives 


*  *  -  '■“j  r.'w, 'jrvy.  vvTdvjir;  J*.Tir;’77,3*.Ti' 


rw  ^  rv* rr.’v.'r.Tv.’' 


Vn  =  Vn_l{l  ♦  r t <eJ6  +  e"J®  -  2)  ♦  r2<eJ*  •*•  e"J*  -  2)}  .  <3.8> 


j  0  _  »  0 

Using  the  identity  e  +e  =2cos©  and  2sin2(0/2)=l-cos0,  Equation 
(3.8)  becomes 


Vn  =  GVn  1  =  {1  -  4r1sin2<2>  -  4r2sin2 ( ^ ) }Vr 


(3.9) 


where  G  is  the  amplification  factor.  Equation  (3.9)  shows  clearly 
that,  if  solutions  are  to  remain  bounded,  we  must  have  IGI<1  for  all  6 
and  4-  This  is  the  stability  criterion  for  the  heat  conduction 


equation. 


For  |G|<1,  we  have 


6  4 

II  -  4r1sin2(2)  “  4r2sin2<->)  <  1 


which  is  true  only  if 


(3.10) 


6  4 

4rtsin2<2>  ♦  4r2sin2<2>  <  2 


for  all  8,  4 


<3. 11) 


The  stability  requirement  is  then 


rt  +  r2  <  2  ’ 


(3.12) 


At  1  1  1 

Rj  (  AC2  +  An2  *  -  2 


(3.13) 


Similarly,  the  equation  for  the  stability  criterion  for  the  substrate 
Is 

1 

Q2<ri  ♦  r2  )  <  -  ,  (3.14) 


or 


Q2At  1  1  1 

"T"  (  AC2  +  An2  *  -  2  * 


(3.15) 


If  Q2  is  less  than  1,  Equation  (3.13)  Is  the  stability  criterion; 
otherwise  Equation  (3. IS)  Is  the  stability  criterion  for  the  current 
problem. 

Based  on  previous  results  in  references  [28,29,30,31,33,34,35],  we 
know  that  high  temperature  and  high  thermal  stresses  occur  In  the 
region  near  the  asperity.  Therefore,  in  that  region  and  in  the  region 
near  the  cavity,  a  very  fine  mesh  must  be  used  to  calculate  accurate 
solutions.  In  the  regions  far  away  from  the  asperity  and  the  cavity,  a 
relatively  coarse  mesh  can  be  used  to  save  computing  time.  This  non- 
uniform  mesh  can  be  transformed  to  a  uniform  mesh  by  using  the  general 
coordinate  transformation  proposed  in  references  [54,55,561.  The  non- 
uniform  mesh  and  general  coordinate  transformation  are  discussed  in 
Appendix  II. 

The  heat  conduction  equation  (3.1)  and  (3.2)  in  the  transformed 
plane  (E,q>  can  be  written  as 


(A.*8.  -  2A2*L  ♦  A,*8.  ♦  A.*!  ♦  A.*8  )/ J2  -  R  f8  ,  3=1,11  (3.16) 

x  tt  2  Cn  nn  *  E  5  n  *  T 


where  J-C_n-  “  £-n-  i*  the  Jacobian  of  the  transformation,  the 
S  n  n  t 


subcripts  <C,q,t,n,T>  denote  partial  derivatives  in  those  coordinates 


and  time,  respectively,  and 


A.  =  U2  ♦  n.2, 

n  n 


<3.17) 


A-  =  c_s_  ♦  n-n-  . 


C  n  t  n 


(3.18) 


a.,  =  t-2  ♦  n-2  . 


t  S 


<3. 19) 


hH  =  <C-A7  -  n-A6)/J 

n  n 


<3.20) 


Ag  =  <n-A6  -  C  A7)/J 


<3. 21) 


A,  =  AJ —  -  2A,C —  ♦  A~C —  , 

6  1  «  2  Cn  3  nn 


<3. 22) 


A,  =  A,n —  ~  2A,n —  ♦  A,n~ 
7  1  «  2  Cn  3  nn 


(3.23) 


In  the  coating  layer,  the  transformed  heat  conduction  equation 


(3.16)  in  the  difference  form  is 


T  I 

♦  <i,j,n>  =  ♦*<i,j,n-l)  ♦  —  AA  , 

Ki 


n=l  <3. 24a) 


,  It  T  2At 

♦  <l.J,n>  =  -  {-♦1<i, j,n-2)  ♦  ^(i.J.n-l)  ♦  -r-  AA)  ,  n>l  <3.24b> 
3  Kx 


where 


v:av; 


AA  =  (<At/A^2)  I4X<i-l, J,n-1>  -  24X<i,J,n-l>  ♦  4X< i*l , J,n-1  >3  - 
-  <A1/2ACAn>C41(i+l, J+l,n-l)-*X<i+l, J-l,n-l)-41<i-l, J+l,n-l>  + 

♦  ♦X(i-1, J-l.n-l)]  ♦  (A3/An2>  t*X(i, j-l,n-l)  -  24I(i,j,n-l>  + 

♦  ♦I(i, J+l,n-l ) ]  ♦  <A4/2AC>  t*X(i+l,3,n-l)  -  ♦X(i-1, J,n-1>3  ♦ 

♦  <A3/2An>  C4X<i,J+l,n-l>  -  #x < 1 . J-l,n-l>]}/J2  .  (3.2S) 


m 


♦*n(i,J+i,n-l)]  ♦  (A4t/2AU  [♦1I(i+l, j,n-i>  -  ♦II(i-l, J,n-i)l+ 

♦<A5/2An)  £*IX(i, j+l.n-1)  -  #II(i, J-l,n-l)3}/J2  .  (3.27) 

At  the  outer  boundaries  of  the  rectangular  region  (excluding  the 
surface),  ♦^<i,J,n)=0  is  the  nominal  value.  The  remaining  conditions, 
on  the  surface,  the  cavity  boundaries  and  at  the  interface,  will  be 
Incorporated  with  an  energy  balance  scheme. 

3.2  Energy  Balance 

The  cavity  boundaries,  the  moving  asperity  and  the  interface  of 
the  medium  are  taken  care  of  with  the  use  of  the  energy  balance 
method  [571. 

(1)  Energy  balance  at  the  interface  (see  Figure  3.2) 

For  material  I  (coating  layer >,  the  heat  fluxes  toward  the  central 
point  P  of  the  element  at  the  interface  from  material  points  W,  R  and  S 
in  the  coating  layer  arc 

T< i-1 , J )  -  T(i,J> 

<W  =  k^Ay/2)  - — -  ,  (3.28) 

T(i+l,J>  -  T(i,J> 

Qr+p  =  k^Ay/2)  - — -  ,  (3.29) 

Ax,+Ax9  T< i , J-l )  -  T(i,J) 

<W  =  v  £  * - £y - •  <3-30> 

where  Q  is  the  heat  flux,  Indexed  by  the  direction. 

For  material  II  (substrate),  the  heat  fluxes  toward  the  point  P 
from  material  points  V,  R  and  N  in  the  substrate  region  are 
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<3. 31  > 


EB 


5 

Qw*p  = 

& 

■> 

* 

Or-»p  = 

■> 

& 

Os->p  ~ 

i 


re 


Sv* 


T(i-l.j)  -  T(i,j) 


Ax, 


T(  i*l , J )  -  T<i,J> 


Ax, 


(3.32) 


II  2 


Ay 


(3.33) 


The  total  heat  flux  going  to  the  interface  point  P(i,j>  Is 


T(i-l,J>  -  T(i,J>  T(i+1, J)  -  T(i,  J> 
Q.u*  *  <k,  ♦  kn)(Ay/2H - - -  ♦  - - - —3  ♦ 


Ax, 


Ax., 


Ax+Ax,  T< 1 , J-l )-T<l, j  >  Ax  +Ax_  T( 1 , J+l )-T( 1 , 1 > 

♦  k  <"V  *> - 7- - ♦  k ln-£> - - -  .( 


I  2 


Ay 


II  2 


Ay 


3.34) 


The  rate  of  change  of  Internal  energy  U  In  the  time  Interval  At  at  the 
point  P<i,J)  Is 


Up  =  ulP  ♦  u2P 


UP  =  <PICI+PIXCII,l(Axl+Ax2>_i  3 


Ay  T(l, J,n>-T(i, J,n-1> 


At 


,  n=l  (3.35a) 


Ur  = 


Ay  3T ( 1 , j , n ) -4T ( 1 , J , n-1 ) +T ( 1 , J , n-2 ) 

<plCl+pliCn)nAxi+Ax2)'4  3  - 


2At 


n>l  (3.35b) 


Conservation  of  energy  requires  that  the  algebraic  sum  of  the  heat 
flowing  into  the  point  P  is  equal  to  the  rate  of  change  of  internal 
energy  at  the  same  point  <QBuaBUp>.  From  conservation  of  energy,  one 


-30- 


vs*: 


J 


can  obtain  the  equation  for  the  continuity  condition  at  the  interface 


4 


(3.39) 


-*<i,j,n-l>]  ♦  I*<i.J+l.n-l>  “  *(i.J,n-l>J 


where  i»fc  =  k^/kj  • 


Details  of  the  energy  balance  method  on  the  other  boundary 
conditions  are  given  in  Appendix  III.  The  dimensionless  form  of  these 


boundary  conditions  are  listed  below. 

<ii>  Energy  balance  on  the  cavity  boundaries  <see  Figure  3.3) 


On  face  AB: 


♦  X<i. J.n)  =  ♦^i.J.n-l)  ♦  AA3  , 


n=l  <3. 40a) 


♦1<i. J.n)  =  \  [-♦1<i, J.n-2)  ♦  4*1(i,J,n-l)  ♦  2AA31  ,  n>l  <3.40b> 


where 


«  —  {[♦I(i-1, j,n-l)  -  24I(i, J.n-1)  ♦  ♦I(i+1, j,n-l)3/A^2  ♦ 


'3  R 


♦  2C*I(i, j-l,n-l>  -  ♦I(i, J.n-1)3/An2) 


(3.41) 


On  face  AC: 


♦  U<i,J,n>  =  ♦II(i,J,n-l>  ♦  AA4  , 


n=l  (3.42a) 


♦”(1,3. n)  =  \  [-♦II(i,J,n-2)+4*II(i,j,n-i)+2AA4)  .  n>l  (3.42b) 


where 


I*. 


ASPERITY 


Figure  3.3  Energy  belence  on  the  cavity  boundary. 


**■»*■««  ■  v  m  vm  i'v 


^  =  {  Ap  t*n<i-l.j.n-l>  -  ♦II<i,j,n-l>]  ♦ 


1  .11 ,  ...  *n . 


An12+AnlAn2 


C*  (i.j-l,n-l>  -  ♦  (i,j.n-l>3  ♦ 


1  .11 


[♦II(i, j+l,n-l)  -  ♦II(i, J.n-1)3} 


An22+AniAn2 


(3.43) 


On  face  BD: 


♦IX(i, J,n)  =  *“(1,3,11-1)  +  AAS  , 


ill 


n=l  (3.44a> 


♦T1(i,j,n)  =  \  C-*II<l.J,n-2>+4*1I(i.J,n-l)+2AA5)  ,  n>l  (3.44b) 


where 


W5  =  ^  (  ^5  [*II(i*l,J,n-l>  -  ♦II(i,j,n-l>3  + 


1  .11,.  -ill- 


C*  <i.j-l,n-l>  -  ♦  <i,J,n-l)J  ♦ 


♦  A  v  1  T—  C*II<1, J+l,n-l)  -  *II(i,J,n-l)]} 
Anj^Ar^Anz 


(3.45) 


On  face  CD: 


*II(l,J,n)  =  *“(l,J,n-l)  ♦  AAe 


n=l  <3. 46a) 


•II<1, J,n)  = 


C-*1I(l,J,n-2)+4*1I(i,J,n-l)+2AA63  ,  n>l  <3. 46b) 
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MV 


kVjOWXi 


wher 


L U  T  T  T  T  T  T 

AA6  =  ^ —  {[*“(i-l, j,n-l)  -  2*1  (i,j,n-l)  «•  ♦“(i+l.j.n-Dl/A*;2  ♦ 
kii 

+  2C*11(i, j*l,n-l>  -  ♦“(i, j.n-l>1}  .  (3.47) 

(111)  Energy  balance  at  the  comer  of  the  cavity  (Figure  3.3) 

The  points  at  the  four  corners  of  the  cavity  are  singularities, 
because  at  each  of  those  four  points  there  are  two  boundary 
conditions,  3T/3x1  =  3T/3x2  =  0,  with  only  one  unknown  T.  However,  by 
applying  an  energy  balance  scheme,  one  can  resolve  such  problems  at 
the  corners.  The  dimensionless  form  for  the  corner  points  are: 

Comer  A: 


♦  <i,j,n>  =  ♦  (  i , J , n-1  >  +  AA.,  , 

1 

*<i,J,n>  =  -  [-♦(!, J,n-2>  +  4*(i,j,n-l>  +  2AA7)  , 


n=l  (3.48a) 


n>l  (3.48h) 


where 


AA .  = 


At 


'  R/2  ♦  *Rtt/4  1  2 

I  k  II 


1+Sk 

[  ~~  [♦(!-!, J, n-1)  -  *( i , j,n-l)]/AS2  + 


+  -  [♦(  i+1 ,  .1  ,n-l  >  -  *(i,  j,n-l)l/A£;2  ♦  ^t*(  i,  ,n-l )  - 


-  *(i,J, n-1) 3/An2  ♦  !♦( i , j-i , n-1 )  -  ♦< i . J,n-1>]/An2)  •  (3.49) 


Corner  B: 


♦<i,j,n>  =  ♦<i,j,n-l)  ♦  AA8  , 


1 

♦  (i, j,n>  =  ^  t  — ♦ < i . j , n-2 )  +  4*<i,j,n-l>  ♦  2AA81  , 
where 


n>l 


AAd 


At 


V2  *  Vii/4 


1 
f 2 


-  #(i,  j,n-l)]/AC2  + 


1+fl.  vu 

♦  -r-*  £♦(!♦!, j,n-l)  -  *<i, j,n-l>]/AC2  +  ^  t*(i,J+l,r 


-  ♦< 1 , j , n-1 ) 1/An2  ♦  C*(i, j-l,n-l>  -  *(i, J,n-1)]/An2} 


Comer  C: 


j,n>  «  *II(i,j,n-l)  +  AA9  ,  n=l 

j,n)  =  \  [-♦II(i,J.n-2>+4*1I<i,J,n-l>+2AA9l  ,  n>l 


where 


AAe 


=  \  p-  |C2#x,<i-l. J.n-1)  -  3*n(i, j,n-l)  ♦  *II<i+l,j 
3  Kli 


(3.50a> 


(3.50b) 


i-l)  - 


(3. SI) 


(3.52a) 


(3.52b) 


,  n-1) 3 /AC 2  + 


♦  [♦II(i,  J-l.n-1)  -  3*”(i,  J,n-1>  ♦  2*n(i,  J*lfn-1> J/An2} .  (3. S3) 


Corner  D: 


♦II(i, J,n)  =  ♦II( i , j  ,n-l )  +  AA10  ,  n=l  <3.S4a> 

j,n)  =  ^  j,n~2)+4iII<i, j,n-l)+2AA103  .  n>l  <3.S4b> 

where 

2 

AA1C  Sr-  { ti11  ( i-1, J,n-1 )  -  3*n(i,  j,n-l>  ♦  2*n<i+l,  j,n-l)  )/AC2+ 
3  Rii 

♦  C*1 1  <  i ,  J-i ,  n-1  > -3*n  <  i ,  J ,  n-1  >+2*X  1  <  i ,  J+l ,  n-1 )  1/An2 }  .  <  3 .  SS ) 

<lv>  Energy  balance  on  the  surface  boundary  (the  moving  asperity) 

The  dimensionless  form  for  the  surface  boundary  condition  is: 

i  i  ^ 

♦  < i,l,n)  =  ♦  <i,l,n*l)  +  <l+h'>— ^  +  AAlt  ,  n=l  <3.56a> 

♦T (i,l,n>  =  \  [-♦I<i,l,n-2)  ♦  4*I(i,l,n-l>  ♦  (l+h’)~^-  ♦  2AA..1  , 

3  KjAn  11 

n=2  (3. S6b) 

where  h'=h/(Ax/2>,  h  is  defined  in  Appendix  III,  and 

AAj t  =  jp  {[♦1<i-l,l,n-l)  -  2*1<i,l,n-l>  «•  ♦* ( i+1 , 1 , n-i ) J/AC2  ♦ 

X 

♦  2t*I(i.2,n-i>  -  ♦I(itl,n-l))/An2J  •  (3.S7) 

Equations  (3.24,  3.26,  3.38,  3.40,  3.42,  3.44,  3.46,  3.48,  3.50, 


** . 


3.S2,  3.54,  3.S6)  constitute  the  general  formulation  of  the  problem 
with  a  complete  set  of  difference  equations  for  the  solutions  of  the 
discrete  temperature  field  {♦<i,j,n>}  at  some  specific  time.  The 
computer  programs  which  are  used  to  compute  the  temperature  field 
solutions  are  given  in  Appendix  IV. 

3.3  Mferical  Results 

Numerical  results  are  obtained  by  using  the  non-uniform  rectangular 
mesh  corresponding  to  different  cases  of  material  properties  and 
geometry.  For  the  surface  layer  of  silicon  carbide,  k  =1.047 
J/cm.°C.s,  *^=0.49  cm2/s,  and  cj=712  J/kg.°C.  For  the  substrate  of 
aluminum,  k  =2.02  J/cm.°C.s,  k  =0.961  cm2/s,  and  c  =917  J/kg.°C.  The 

II  II  II 

other  numerical  parameters  on  the  asperity  and  the  cavity  are:  v=15 

m/s,  w=10a,  H=1.2a,  b=1.9a.  d=0.3a,  e=0.5a,  a=lmm,  the  smallest  A£  and 

An  are  0.02  and  0.01  respectively,  and  At=0.01.  In  the  limiting  case 

of  no  cavity,  the  maximum  dimensionless  temperature  at  the  surface  ol 

the  coated  media  was  found  to  be  0.124  by  using  the  Fourier  transform 

method  [33,  34],  The  result  at  the  same  point  by  the  current  finite 

difference  formulation  is  0.123.  The  error  is  less  than  12.  The 

numerical  scheme  is  therefore  confirmed  by  the  benchmark  problem. 

The  solutions  for  a  single  material  with  and  without  a  cavity  would 
then  be  compared  with  two  limiting  cases.  For  the  first  case,  the 
cavity  is  located  entirely  in  the  surface  layer,  Figure  3.4a.  In  the 
second  case,  the  top  edge  of  the  cavity  is  at  the  layer/substrate 
interface,  Figure  3.4b.  The  solutions  for  the  single  material  without 
and  with  a  cavity  are  designated  as  the  third  and  the  fourth  cases, 
respectively,  included  for  the  purpose  of  comparison.  Different  cases 
of  the  temperature  field  solutions  are  given  in  Table  1.  In  Case  1, 
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Table  1 


I 


e* 


case 

kI 

kn 

<pc>r 

<pc’n 

KI 

K 

II 

L 

1A 

lkl 

lkn 

1 < pc )  1 

1,pc,„ 

1KI 

lKII 

0.04 

IB 

lkI 

lkn 

1  <  pc  )  z 

1<pc>II 

1KI 

lKII 

0.06 

1C 

lkI 

lkn 

Kpclji 

2ki 

1KII 

0.04 

ID 

K 

M 

M 

A 

T-4 

4<pc,r 

1<pc>ir 

1ki 

lKII 

0.04 

IE 

lkl 

lkn 

1 <  pc )  x 

1<PC,H 

lKi 

lKII 

0.1 

2A 

lki 

lkn 

KpOj 

1<pc,u 

lKi 

lKII 

0.04 

2B 

lkl 

2kr 

l<pc)j 

2(pc>i 

lKx 

lKI 

0.04 

2C 

lkI 

2ki 

l(pc>i 

j<pp>, 

lKr 

lKI 

0.04 

2D 

lkI 

lki 

1 <pc)j 

|<ppi, 

lKi 

2KI 

0.04 

2E 

lki 

2k: 

l(pc)i 

1 (pc)j 

lKi 

2ki 

0.04 

3 

lki 

Kpclj 

lKi 

4 

lk 

l(pc) 

lK 

0.04 

cavity  location 
coating  layer 
coating  layer 
coating  layer 
coating  layer 
coating  layer 
interface 
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•  Base  materials  for  the  coating  layer  and  the  substrate  are  silicon 
carbide  and  aluminum,  respectively. 


three  different  values  of  the  ligament  thickness:  0.04  (Case  1A>,  0.06 
(Case  IB),  and  0.1  (Case  IE),  are  used  to  illustrate  the  effect  of  the 
ligament  volume.  The  temperature  fields  at  two  depths,  for  all  three 
Cases  1A,  IB  and  IE  are  shown  in  Figure  3.S.  In  the  figure,  the  cavity 
width  is  from  E=1.6  to  2.2.  The  asperity  position,  at  the  dimension¬ 
less  time  x=1.04,  respresenting  the  worst  case,  is  from  6=1.2  to  2.2. 
The  relative  positions  of  the  asperity  and  the  cavity  is  shown  in 
Figure  3.6. 

The  Case  1A  then  is  compared  with  Case  2  of  the  same  ligament 
thickness  for  which  the  top  edge  of  the  cavity  is  at  the  layer/ 
substrate  interface.  The  effect  of  the  relative  position  of  the 
cavity  to  the  interface  is  shown  in  Figure  3.7.  It  is  noticed  that 
when  the  top  edge  of  the  cavity  is  at  the  interface,  the  temperature 
field  in  the  region  immediately  on  the  trailing  edge  of  the  asperity 
will  be  affected  by  the  substrate  material. 

The  effect  of  the  heat  capacity  and  thermal  conductivity  of  the 

surface  layer  for  Case  1A  is  shown  in  Figure  3.8.  The  figure  shows 

the  original  value  as  Case  1A.  Case  1C  represents  a  reduction  of 

thermal  capacity  of  the  surface  layer  by  half.  Case  ID  shows  the 

result  of  an  increase  in  thermal  conductivity  of  the  surface  layer  by 

7SZ.  The  thermal  conductivity  of  the  surface  layer  is  shown  to  have 

little  effect  on  the  nondimens ional  surface  temperature.  But  the  real 

temperature  field,  T=q  a$/k  t  Is  lowered  with  Increasing  thermal 

o  i 

conductivity  k  . 

Figures  3.9  and  3.10  illustrate  the  effect  of  a  cavity  on  the 
direction  of  heat  flux.  The  figures  show  the  nondimensional  heat  flux 
components  in  6  and  n  directions  of  a  single  material  without  a  cavity 
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Figure  3.7  Dimensionless  temperature  field  (cases  1A,2A,3> 


Figure  3.10  Dlaenslonless  heat  flux  in  n  direction  (cases  1A,3> 
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(Case  3)  and  a  layered  medium  with  a  cavity  (Case  1A>.  From  the 
figures,  it  is  observed  that,  with  no  cavity,  the  heat  flux  at  £=2.2, 
and  n=0.04  has  a  magnitude  of  0.7  at  an  angle  of  82°  to  the  wear 
surface.  With  a  cavity,  at  the  same  location,  the  magnitude  is 
Increased  to  l.S  at  an  angle  23°  to  the  wear  surface.  Hence,  the 
existence  of  the  cavity  will  increase  the  heat  flux  tremendously, 
especially  in  the  £  direction  near  the  upper  trailing  corner  of  the 
cavity.  Figures  3.9  and  3.10  demonstrate  not  only  an  Increase  in 
magnitude  of  the  heat  flux,  hence  the  temperature  gradient,  but  also 
the  flux  at  a  more  oblique  angle  to  the  wear  surface. 

Ju  1281  has  studied  the  effect  of  thermal  properties  of  a  single 
material  subjected  to  the  high-speed  asperity  excitation.  It  was 
pointed  out  that  thermal  conductivity  (k>  and  thermal  capacity  (pc)  are 
the  parameters  controlling  the  temperature  field.  For  layered  media,  a 
similar  effect  was  found  by  Ju  and  Liu  1351.  For  the  case  of  a  layered 
medium  with  a  cavity,  the  thermal  property  variation  in  the  coating 
layer  can  be  accordingly  extrapolated.  It  is  the  effect  of  the 
substrate  in  the  neighborhood  of  the  cavity  that  would  be  influential 
in  determining  the  temperature  field  in  the  critical  region.  The 
effect  of  thermal  property  variation  for  the  substrate  is  therefore 
studied  numerically  for  the  Case  2,  for  which  the  coating/ substrate 
interface  is  at  the  top  edge  of  the  cavity.  For  this  case,  the  thermal 
properties  of  the  substrate  will  be  of  immediate  influence  to  the 
temperature  field  in  the  vicinity  of  the  top  trailing  corner  of  the 
cavity.  For  the  purpose  of  demonstrating  the  individual  effect,  a 
benchmark  case  is  chosen  for  comparison  in  which  both  the  coating  and 
the  substrate  are  of  silicon  carbide,  <kj=kIi=1.047  J/cm.°C.s, 
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ki=kii=0,49  c“2/s’  pici=piicii=2137  J/c®3-C>t  designated  as  Case  4. 
Figure  3.11  illustrates  the  temperature  field  near  the  top  surface  and 
at  the  coating/substrate  interface  for  cases  with  marked  changes  in 
thermal  properties  from  those  given  in  Case  4.  Case  2B  shows  no  change 
in  the  substrate  diffusivity,  but  both  the  thermal  conductivity  and  the 
thermal  capacity  are  doubled.  The  ensuing  improved  conductivity  and 
capacity  in  the  substrate  allow  a  significant  heat  flow  into  the 
substrate;  thus  a  high  temperature  gradient  is  also  found.  The  Case 
2C,  at  the  same  diffusivity,  but  with  both  the  thermal  conductivity  and 
the  thermal  capacity  halved,  shows  a  reduced  heat  flow  into  the 
substrate,  with  a  corresponding  low  temperature  gradient.  Cases  2D  and 
2E,  with  doubled  diffusivity,  but  with  half  capacity  and  double 
conductivity,  respectively,  showed  reduced  and  increased  heat  flow  into 
the  substrate,  respectively.  The  heat  flux,  being  proportional  to  the 
temperature  gradient,  is  illustrated  in  Figures  3.12  and  3.13  for  the 
surface  region  and  at  the  interface. 

Figure  3.14  shows  the  transient  temperature  for  Case  1A  (cavity  in 
the  coating  and  ligament  thickness  of  0.04>  in  comparison  to  the  case 
of  a  single  material  without  cavity  (Case  3>.  The  dimensionless 
temperature,  ♦=Tkj/qoa,  plotted  against  dimensionless  time,  *=vt/a,  at 
the  surface  and  at  the  ligament  depth,  rp0.04,  for  the  position  K-2.2, 
where  the  temperature  is  maximum  in  the  vicinity  of  the  cavity.  It  is 
shown  that,  before  the  asperity  reaches  the  point,  the  temperature  is 
low.  Then  the  surface  temperature  increases  and  reaches  a  maximum 
when  the  asperity  just  passes  over  the  trailing  edge  of  the  cavity. 
After  the  passing  of  the  asperity,  the  temperature  at  the  trailing 
corner  of  the  cavity  drops  again. 
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Figure  3.14  Dluenslonless  te^jeratwe  field  (cases  1A,3). 


CHAPTER  4 


STRESS  SOLUTION 

Ju  et  al  128,29,30,31,32,33,34]  established  that,  for  a  moderately 
high-speed  asperity  excitation,  the  thermal  stress  effect  dominates  the 
stress  field  and  eventually  leads  to  failure  in  the  no-cavity  case. 

Liu  CS8]  in  his  thesis  also  showed  that,  if  the  asperity  speed  is 
larger  than  0.127  m/s  <S  in/s),  the  thermal  stress  dominates  the 
failure,  and  the  mechanical  stress  becomes  less  Important.  However, 
the  mechanical  stress  may  not  be  trivial  when  a  cavity  exists. 
Therefore,  both  the  mechanical  and  thermal  stress  field  will  be 
presented  in  this  chapter. 

The  thermoelastic  Navier's  equations  and  the  Hooke's  law  in 
dimensionless  form  are: 
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where  u$  and  v^  are  respectively  the  dimensionless  stress  modulated  by 
the  average  pressure  P0;  M=(v/c2>  is  the  Mach  number;  V  is  the 
asperity  speed;  +  2y3)/pnc2i  N^X^/p^c,,;  Ng^/P^c.,; 

bI=<3V2p3>/pxi;  V(q°aV/ki;  6=Vpii’  ci  =  t(Xii+  2pn)/pii]l/2’ 

c2  -  (pii/pi t ) 1^2  are  the  dilatational  and  shear  wave  speed, 

respectively;  L,  q,  -t,  X,  p„,  a.,  p_  and  $3  were  defined  in  Chapters  2 

P  P  P 

and  3. 

Equations  for  the  mechanical  stress  field  are  the  same  as  the 
thermal  stress  field  except  that  there  is  no  temperature  effect. 
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4.1  Perturbation  Method 

For  hard  wear  materials,  such  as  Stellite  III,  the  Mach  number  M  is 
of  the  order  of  10_3.  Since  M2  is  a  parameter  which  is  sufficiently 
small,  Equations  (4.1)  and  <4.2>  can  be  solved  by  the  perturbation 
method. 

Let  the  solutions  to  (4.1)  and  (4.2)  be  expressed  as  a  power 
series  of  e=M2;  that  is 

u£(E ,  q,T,e  >  =  u§<E,q,T>  +  euf<E,q,T:>  +  e2u§(E,q,T)  +  •••  .  (4.6) 

v£(t;  ,q,T,e )  -  v^(L,q,T)  +  ev^(E,q,t>  +  e2v^<E,q,'t)  +  ...,  (4.7) 

when  equations  (4.6)  and  (4.7)  are  substituted  into  equations  (4.1) 
and  (4.2),  the  terms  with  the  same  power  of  e  are  grouped,  leading  to 
a  set  of  equations  for  u§,  u^,  u^,  ...,  and  v^,  v^,  v^,  ...,  as 
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follows : 


For  the  terms  of  the  zeroth  order  in  e, 
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For  the  terms  of  the  first  order  in  e, 


3u§ 


—  <N,~rr)  + 


35 


35 


3v§  3  3v§  3  3u§ 

(N0-~)  +  —  <N«-~>  +  rr  (N..-rr>  = 


35  2  3n  3n  3  35  9n 


3q 


=  6 


3x2  ’ 


(4.10) 


3  dv§  3  3u§  3  3u£  3  3v? 

35  <N3_35>  +  3rj  <Nz  35 J  +  35  <N3_an>  +  a^  <Ni  an> 


3q  3q 


3n 


=  6 


3t2  * 


(4.11) 


For  the  higher  order  solutions,  it  is  evident  that  the  equations  are 
recursive.  Accordingly,  the  recurrence  formulas  can  be  written  as: 
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where  6i0  is  Kronecker  delta. 


Similarly,  we  can  obtain  a  set  of  equations  for  the  stress  field  as 


follows : 
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The  recurrence  formulas  for  the  stress  field  are 
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The  boundary  conditions  are  as  follows: 


For  the  zeroth  order  solutions, 


<i)  the  surface  (c(t)/a  <  £  <Cc(t)/a  +  11,  n=0>  is  traction  prescribed 
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where 


_  /  P’<£>  ,  for  the  mechanical  stress  field 
V  0  .  for  the  thermal  stress  field 


(ii)  the  regularity  conditions  at  infinity,  are 


ug,  vg,  °o|r|,  =  0 


(iii>  the  continuity  conditions,  at  n=0,  are 


11I  =  |]I  ^ 
uO  U0  1 


yl  =  vjl 
vo  V0  * 


°o *  =  °alI  , 

Sn  Sn 


°0*  =  OqII 

nn  nn 


i’l  l*!,* 


(4.23) 


(4.24) 


(4.25) 


(4.26) 


(4.27) 


(iv)  the  cavity  boundary  conditions  are  traction  free. 


For  the  nth  order  solutions,  the  surface  is  traction  free  and  the 
remaining  boundary  conditions  are  the  same  as  for  the  zeroth  order 
solutions.  The  solutions  of  each  perturbative  order  can  be  obtained 
by  applying  the  finite  difference  method. 

4.2  Difference  Formulation 

Because  of  the  complexity  of  the  geometry  and  the  boundary 
conditions,  the  finite  difference  method  is  considered  more  appropriate 
than  the  transform  method,  which  was  used  in  the  cases  without  a 
cavity.  In  this  section,  only  the  zeroth  order  solutions  of  the 
thermal  stress  field  will  be  discussed  in  detail,  the  solutions  of 
higher  order  and  the  mechanical  stress  field  can  be  obtained  similarly. 

In  the  finite  difference  method,  the  semi-infinite  body  is  replaced 
by  a  sufficiently  large  rectangular  region  (Figure  4.1),  and  a  non- 
uniform  mesh  must  be  used  as  we  stated  in  Chapter  3.  The  non-uniform 
mesh  is  transformed  to  the  uniform  mesh  by  applying  the  general 
coordinate  transformation  (Appendix  II).  The  stress  field  can  then  be 
solved  in  the  transformed  plane  (computational  plane)  (£,n>-  The 
finite  difference  form  of  the  thermoelastic  Navier's  equations  (4.8) 
and  (4.9)  in  the  computational  plane  (CfH>  can  he  written  as: 
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The  traction  surface  boundary  conditions  are  expressed  in  terms  of 
the  displacements: 

-  vI(i-l,l,n)/(2C-At)  -  3ur( i ,1 ,n)/(2r)-An>  ♦  2uI(l,2,n)/(n-An>  " 
t  n  n 


<4.48) 


-  u^i.a.nJ/CZn-An)  +  v^i+l.l.n)/^-^)  =  o  , 

n  £ 

and 

-  N,<i,l)ur<i-l,l,n)/<2UAt;>  -  3N,  <  i  .  1  >  v1  <  i ,  1 ,  n>/<  2n-An>  + 

£  t  n 

+  2N. < i,l )vI< i ,2,n)/<n-An>  -  N. ( i , 1 )vx < i , 3, n)/ < 2n-An>  + 

n  n 

b?Y  j 

+  N2<  i ,  1  )ur  <  i+1 ,1 ,  n)/<  2E_Ai; )  =  ~L~~  ^(1,1, n>  .  <4. 49) 

S  C2 

The  traction  free  boundary  conditions  on  the  cavity  (Figure  3.3) 

are: 

On  face  AB : 

-  vr(  i-1 ,  j  , n)/(2£_A£  >  +  uI(i,j-2,n)/<  2n-An )  -  Zu^i,  j-l,n)/<n-An> 

K  n  n 

+  3uI(i, j,n)/(2o-An>  ♦  v1 < i+1 , j , n)/( 2^_A^ >  =  0  ,  <4. SO) 

n  S 

and 

-  Nz(i,  j>uI(i-l,  j,n)/<2{;_AF)  +  N,  <  i ,  J  >v*<  i ,  j-2,  n)/<  2n-An>  - 

-  2N.(i,  j)vI(i,  j-l,n)/<n-An>  +  3N.<i,J)vI(i,  j ,  n)/(  2r)-Ar)  >  + 

n  n 

+  N2(i,  J)uI<i+l,  j,n)/<2£;-AF)  =  ♦1<i,j,n)  , 


(4.51) 


On  face  AC: 


-  N,(i,  j)uII<i-l,  j,n)/<{._AO  -  N2<  i ,  j  >vx  1  <  i ,  j-1 ,  n>/(2n-An>  + 
1  K  n 


+  Nl<i,j)uII(i,j,n)/<  S-A£ )  +  N2( i . j ) v11 ( i , j+1 , n>/( 2n-An> 


biiYn  .  , 

___  ♦II(i, J,n>  , 


(4.52) 


-  vII<i-l,  i,n)/(f;_^)  -  uIX<i,  j-l,n>/(2n-An)  +  v11^,  j,n)/(f;-At)  ♦ 

c  n  t 


+  uir<i, j+l,n)/(2r)_An)  =  0  , 

n 


(4.53) 


On  face  BD: 


-  N2<i, j)vII<i, j-l,n)/(2n-An>  -  Nj ( i , j >uIX< i , j , n>/U.-AC >  + 


+  N2(i,  J  >vII(i,  j+l,n)/(2ri_An>  +  Nj  ( i ,  j  )uIX  ( i+1 ,  j , n)/<t-A^>  = 


b?irll  TT  .  , 

c2 


(4.54) 


-  urr(i, j-l,n)/(2n-An>  “  vir<i, j,n)/<C_AC)  + 
H  * 


♦  u11  (  i ,  1+1 ,  n)/(2n-An>  +  vII(i+l,  j,n)/<t;_AC)  =  0 


(4.SS) 


'Cmr'l  WT  ' 


On  face  CD: 


-  vII(i-l, j,n)/(2E_AC)  -  3u* 1  < i , j , n)/( 2n-Ar| >  + 

t  n 


+  2U11 < i, j+1 , n)/< n-Aq >  -  u11 ( i , j+2, n)/(2p_Ap>  + 

n  n 


+  v1 1  ( i+1 , j , n>/( 2£_A£ )  =  0  , 


(4.56) 


-  N2(l, j)uIX(i-l, j,n)/(2UA£>  -  3N. < i , j )vri ( i , j , n>/( 2p_Ap>  + 

5  n 


+  2N.< i , j )vrr < i , j+l,n)/(p_Ap)  -  N, ( i , j )vl 1 ( i , j+1 , n)/(2p_Ap )  + 

n  n 


b?-iYr  t 

+  N2<l(j>uzz<l+l,J,n>/<2E-/£>  =  ~jr^k  ♦II<ilj,n>, 

S  C2 


(4.57) 


where  $P(i,j,n>,  3$P(i,  j,n)/3t,  and  3$^( i , j , n)/3p  are  input  data 


obtained  from  the  temperature  field  solutions  (re:  Chapter  3). 


4.3  Cavity  Corners  Singularities 


When  values  of  a  solution  of  a  boundary-value  problem  or  its 


derivatives  approach  infinity  at  points,  lines,  or  surfaces  in  the 


domain,  the  solution  is  said  to  possess  singularities  at  these  places. 


The  approximation  of  functions  with  singularities  presents  some  serious 


numerical  difficulties.  Nevertheless,  calculation  of  solutions  with 


singularities  is  extremely  important;  such  problems  arise  in  fracture 


mechanics,  various  flow  phenomena,  heat  conduction  problems,  and  in 


fact,  in  any  boundary-value  problem  in  which  strong  irregularities 
occur  in  one  or  more  of  the  following:  (a)  the  geometry  of  the  domain, 
(b)  the  coefficients  in  the  governing  differential  equation,  or  <c)  the 
prescribed  functions,  and  so  on. 

Despite  the  difficulties ,  numerical  methods  can  be  devised  that 
yield  excellent  results  for  singular  problems.  Basically,  there  are 
two  general  ways  the  problem  can  be  approached: 

Nonuniform  Meshes:  This  means  that  a  finer  gradation  of  the  mesh 
is  used  in  the  neighborhood  of  singular  points  in  order  to  capture 
large  changes  in  the  gradients  of  the  solution  nearby.  This  is  often 
a  straight  forward  and  effective  way  to  handle  singularities  and  it 
requires  no  special  modification  of  the  code  or  special  elements,  but 
it  may  be  expensive  owing  to  the  necessity  of  a  large  number  of  grid 
points. 

Special  Singular  Elements:  The  scheme  in  this  case  is  to  devise 
special  elements  in  which  the  approximation  simulates  the  diverging 
rate  in  elements  in  the  vicinity  of  the  singular  point.  However,  this 
method  can  be  used  only  when  the  behavior  of  the  singularity  is  known. 
The  procedure  of  this  method  is  to  assume  a  series  which  consists  of 
both  the  regular  terms  and  the  singular  terms.  For  thermomechanical 
problems,  the  series  form  of  the  asymptotic  expansion  can  be  written  in 
the  form  r59,60,61 ,62,63] : 

u(r,0)  =  regular  term  +  E  AnrnT|/^  f(0)  ,  (4.S8) 

where  r,  0,  and  C  are  defined  in  Figure  4.2.  Indeed,  when  n/C  is  not 
an  integer,  the  derivatives  of  the  leading  term  in  the  singular  part  of 


u  may  become  unbounded.  The  order  of  the  singularity  increases  as  C 
increases.  If  n/Q  <  1,  P  is  referred  to  as  a  reentrant  corner  and  the 
first  derivatives  of  u  are  then  unbounded  as  r-»0.  For  the  present 
problem,  C=3n/2. 

For  the  current  problem,  the  stress  singularity  at  the  cavity 
corner  can  be  resolved  by  using  the  results  of  Williams  [641  and  Sih 
F6S1.  The  series  form  for  the  displacements  in  the  neighborhood  of  the 
cavity  corner  are: 

u(r,0)  =  a1r2/3cos(20/3)  +  a2r2/,3sin<  20/3 )  +  a3r2/3cos(  40/3 )  + 


+  a„r2^3sin(40/3>  +  a-r<*/3cos(40/3)  +  a,r1|^3sln(40/3)  + 


+  regular  term  , 


(4.59) 


v(r,0>  =  b1r2/3cos(20/3)  +  b2r2/3sin( 20/3 )  +  b3r2/3cos(40/3>  + 


+  b„r2/3sin<40/3)  +  b«W3cos<40/3>  ♦  bAr<‘/3cos(40/3)  + 


+  regular  term  . 


(4.60) 


In  the  special  elements,  parts  of  the  coefficients  of  the  series  can 
be  determined  by  substituting  equations  (4.59)  and  (4.60)  into 
Navier's  Equations  (4.8)  and  (4.9)  and  the  traction  free  conditions  on 
the  cavity.  The  remaining  coefficients  can  then  be  solved  by  using 
the  difference  form  (4.28)  and  (4.29).  As  illustrated  in  Figure  4.3 


for  one  specific  corner,  one  can  use  Equations  <4. 8)  and  <4. 9)  at 
points  3  and  7  and  the  traction  free  condition  at  point  5  incorpora¬ 
ting  the  singular  behavior;  then  one  can  use  the  difference  form 
(4.28)  and  (4.29)  at  points  1,  2,  4,  6  to  determine  the  other 
coefficients  of  the  series. 

Equations  (4.28,  4.29,  4.48,  4.49,  4.50,  4.51,  4.52,  4. S3,  4.54, 
4.SS,  4.S6,  4.S7)  and  the  special  elements  compose  a  complete  set  of 
difference  equations  for  finding  the  stress  field.  For  the  zeroth 
order  solutions,  we  have  to  solve  a  set  of  simultaneous  algebraic 
equations,  which  can  be  separated  into  two  groups,  depending  on  whether 
the  coefficient  matrix  is  dense  (few  zero  elements)  or  the  coefficient 
matrix  is  sparse  (many  zero  elements).  The  two  commonly  used  methods 
of  solving  simultaneous  algebraic  equations  include  the  direct  method 
and  the  iterative  method  C66,  67,  68,  69,  70,  711. 

Figure  4.4  shows  the  element  pattern  of  the  matrix  for  the  zeroth 
order  solutions.  It  is  a  largo,  banded,  but  unsymmetric  matrix. 

Because  of  the  dimension  of  the  matrix  <=  4400x4400),  it  is  almost 
impossible  (too  expensive)  to  store  all  of  the  elements.  Fortunately, 
the  matrix  is  banded,  therefore  we  can  only  store  the  elements  inside 
the  bandwidth  by  using  the  one-dimensional  array  as  shown  in  Figure 
4.5,  and  then  using  Gauss  elimination  to  solve  the  system.  The 
computer  programs  for  solving  the  thermal  stress  field  are  given  in 
Appendix  IV. 

4.4  llimerical  Results 

For  the  numerical  examples.  Stellite  III  is  used  as  the  base 
material.  The  mechanical  and  thermal  properties  of  stellite  III  are: 
E=240xl03  Mpa,  v=0.28S,  p=8.3x!03  kg/m3,  K=9.7  J/m-°K-s,  x=2.77xl0-6 


MZ 
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ra2/s,  a=11.3xlO“6  m/m°K,  and  pf=0.5.  For  this  problem,  the  smallest 
AC  and  An  in  the  stress  field  are  0.02  and  0.006,  respectively.  The 
total  grid  points  in  £  and  n  directions  are  67  and  35.  The  other 
important  numerical  parameters  are:  V=15  m/s,  w=30a,  d=0.3a,  e=0.5a, 
and  a=l  mm. 

In  the  numerical  results,  the  effect  of  the  cavity  location  and  the 
effects  of  the  mechanical  and  thermal  properties  on  the  stress  field 
are  also  studied.  All  figures  (Figures  4.7  through  4.24)  are  plotted 
for  the  worst  case  of  the  asperity  position,  that  is  when  the  asperity 
is  right  over  the  cavity  or  when  the  trailing  edges  of  the  asperity  and 
the  cavity  are  aligned  as  shown  in  Figure  4.6  (asperity  position  in 
from  £=-0. 7  to  £=0.3;  cavity  location  is  from  £=“0.3  to  £=0.3). 

When  the  cavity  is  located  entirely  in  the  surface  layer,  because 
the  coating  layer  is  thick,  the  effect  is  similar  to  the  effect  of  a 
single  material  [33,34,3S].  Figures  4.7  to  4.12  plot  the  thermal 
principal  stresses  along  the  asperity  traverse  direction  at  the 
critical  depth  for  the  cases  of  a  single  material  with  a  cavity. 
Different  cases  of  a  single  material  with  a  cavity  are  tabulated  in 
Table  2.  Figure  4.7  compares  the  dimensionless  principal  thermal 
stress  of  the  single  material  with  (case  1A)  and  without  (case  2)  a 
cavity.  The  maximum  dimensionless  tensile  thermal  stress  is  0.98  for 
no-cavity  case  occurring  at  a  depth  q=0.16,  while  it  is  5.9  at  a  depth 
of  0.088  for  the  medium  with  a  cavity  at  a  ligament  thickness  0.094. 

The  location  of  the  cavity  from  the  wear  surface,  as  Indicated  by 
the  ligament  thickness,  influences  the  temperature  field  in  the  total 
volume  available  for  heat  content  generated  by  frictional  heating.  As 
a  consequence,  the  thermal  stress  state  is  strongly  affected.  Figure 
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Table  2 


case 

k 

pc 

E 

a 

L 

cavity 

stress  computed 

1A 

lk 

lpc 

IE 

la 

0.094 

Yes 

thermal  stress 

IB 

lk 

lpc 

IE 

la 

0.06 

Yes 

thermal  stress 

1C 

lk 

lpc 

IE 

la 

0.122 

Yes 

thermal  stress 

ID 

lk 

lpc 

IE 

la 

0.159 

Yes 

thermal  stress 

2 

lk 

lpc 

IE 

la 

No 

thermal  stress 

3A 

!k 

lpc 

IE 

la 

0.094 

Yes 

thermal  stress 

3B 

2k 

lpc 

IE 

la 

0.094 

Yes 

thermal  stress 

4A 

k 

2  pc 

IE 

la 

0.094 

Yes 

thermal  stress 

4B 

k 

1 

2pC 

IE 

la 

0.094 

Yes 

thermal  stress 

SA 

k 

pc 

3E 

la 

0.094 

Yes 

thermal  stress 

SB 

k 

pc 

2E 

la 

0.094 

Yes 

thermal  stress 

SC 

k 

pc 

2E 

la 

0.094 

Yes 

thermal  stress 

6A 

k 

pc 

E 

2a 

0.094 

Yes 

thermal  stress 

6B 

k 

pc 

E 

b 

0.094 

Yes 

thermal  stress 

•  Base  material  is  Stellite  III. 


Figure  4.7  Thermal  principal  stress  (cases  1A,2). 


4.8  shows  the  effect  of  various  cavity  locations  on  the  thermal  stress 
field.  Ju  et  al  [27,32,S8]  showed  that,  without  a  cavity,  the  critical 
depth  at  which  the  tensile  thermal  stress  reaches  a  maximum  is  rpO.16 
for  a  moving  line  asperity  excitation  over  Stellite  III.  However,  with 
a  near-surface  cavity,  both  the  temperature  distribution  and  its 
gradients  are  different,  thus  changing  the  critical  depth.  The  maximum 
tensile  stress  is  optimized  with  respect  to  various  ligament  thickness. 
The  worst  case  for  the  maximum  tensile  stress  for  Stellite  III  is 
obtained  when  the  top  edge  of  the  cavity  is  at  n=0.094. 

Figure  4.9  presents  the  effect  of  the  thermal  conductivity.  Case 
1A  is  of  Stellite  III.  Case  3A  shows  the  effect  that  the  thermal 
conductivity  is  reduced  by  half,  while  case  3B  demonstrates  the  effect 
that  the  thermal  conductivity  is  twice  that  of  Stellite  III.  Figure 

4.9  establishes  that  the  principal  thermal  stress  increases  with 
decreasing  thermal  conductivity.  In  Figure  4.10,  cases  4A  and  4B  show 
the  results  of  doubling  and  halving  the  thermal  capacity,  respectively. 
We  observe  that  decreasing  heat  capacity  will  increase  thermal  stress. 
Figures  4.11  and  4.12  demonstrate  the  effects  of  Young's  modulus  and 
the  coefficient  of  thermal  expansion.  In  Figure  4.11,  Young's  moduli 
for  cases  5A  and  SB  are,  respectively,  three  times  and  one-half  that  of 
Stellite  III.  In  Figure  4.12,  the  thermal  expansion  coefficients  for 
cases  6A  and  6B  are,  respectively,  twice  and  one-half  that  of  Stellite 
III.  These  two  figures  clearly  show  that  increasing  either  Young’s 
modulus  or  the  thermal  expansion  coefficient  induces  higher  thermal 
stress . 

When  the  top  edge  of  the  cavity  is  at  the  interface,  both  the 
coating  layer  and  the  substrate  will  influence  the  stress  field. 


Figure  4.9  Thermal  principal  stress  (cases  1A, 3A, 3B) .  L=0. 094,  and  n=0.06  above 


principal  stress  (cases  1A,4A,4B>.  L=0.094,  and  n=0-06  above 


.094,  and  n-0.06  above 
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Figure  4.12  Thenal  principal  stress  (cases  1A,6A,6B> .  L= 0.094,  and  n=0.06  above 


Figures  4.13  to  4.24  show  the  results  of  the  cases  in  which  the  top 
edge  of  the  cavity  is  at  the  interface.  Different  cases  of  a  layered 
medium  with  a  cavity  are  listed  in  Table  3. 

In  the  case  of  a  medium  with  no-cavity,  the  effect  of  the 
mechanical  stress  field  is  small  enough  to  be  neglected.  When  a  cavity 
exists,  the  effect  of  the  mechanical  stress  field  is  no  longer 
negligible.  Figure  4.13  plots  the  principal  thermal  stress  field  (case 
7A),  mechanical  stress  field  (case  7B),  and  combined  stress  field  (case' 
70.  In  this  figure,  the  material  of  the  substrate  is  Stellite  III, 
and  the  material  properties  of  the  coating  layer  are  the  same  as 
Stellite  III  except  that  Young's  modulus  Is  twice  that  of  Stellite  III. 
This  figure  establishes  that  the  tensile  thermal  stress  is  larger  than 
the  tensile  mechanical  stress.  However,  the  mechanical  stress  field  is 
not  so  small  that  we  can  neglect  it  as  indicated  in  the  no-cavity  case. 

The  effect  of  the  cavity  location  on  the  thermal  stress  field  for  a 
layered  medium  is  presented  in  Figure  4.14.  From  this  figure,  one  can 
see  that  the  maximum  tensile  stress  occurs  when  the  ligament  thickness 
L=0.094,  which  is  the  same  value  as  in  the  case  of  a  single  material 
with  a  cavity.  Figures  4. IS  and  4. IS  present  the  effects  of  Young’s 
modulus  of  the  coating  layer  and  the  substrate,  respectively.  In 
Figure  4. IS,  the  material  of  the  substrate  is  Stellite  III  for  all 
cases.  Young's  modulus  of  the  surface  layer  for  different  cases  is: 
case  1A  is  the  same  as  Stellite  III;  case  7A  is  twice  that  of  Sellite 
III,  case  9A  and  case  9B  are,  respectively,  three  times  and  one-half 
that  of  Stellite  III.  From  this  figure,  one  can  see  that  the  principal 
thermal  stress  field  is  strongly  Influenced  by  the  Young's  modulus  of 
the  coating  layer;  increasing  Young's  modulus  of  the  coating  layer  will 


Table  3 


case 

ki 

kn 

El 

En 

aI 

an 

L 

cavity 

stress  computed 

7A 

Ikj 

Ikn 

2EX 

1EIX 

lax 

laXI 

0.094 

Yes 

thermal 

stress 

7B 

lkj 

lkXI 

2EX 

IE  n 

lOj 

l«i  i 

0.094 

Yes 

mech. 

stress 

7C 

lkx 

lkXI 

2EX 

lEir 

iOj 

laXI 

0.094 

Yes 

combined  stress 

8A 

Ik  i 

lkn 

2EX 

1EIX 

l°i 

lan 

0.06 

Yes 

thermal 

stress 

8B 

lkj 

lku 

2Ej 

lEn 

lax 

laIX 

0.07 

Yes 

thermal 

stress 

8C 

lkt 

lkn 

2EX 

lEn 

lai 

laIX 

0.12 

Yes 

thermal 

stress 

8D 

Iki 

lkn 

2EX 

1Eu 

l°i 

la„ 

0.159 

Yes 

thermal 

stress 

9A 

lkx 

lkjj 

~E 

2E* 

1En 

lax 

1°ii 

0.094 

Yes 

thermal 

stress 

9B 

lk, 

lkn 

3EX 

1E„ 

l«x 

1°ii 

0.094 

Yes 

thermal 

stress 

10A 

lkx 

lkjj 

lEj 

lr 

sEir 

lar 

lan 

0.094 

Yes 

thermal 

stress 

10B 

lkx 

lkn 

lEi 

2Eix 

laj 

laxx 

0.094 

Yes 

thermal 

stress 

IOC 

Ik  i 

lkn 

IE  I 

SEn 

lax 

1°ii 

0.094 

Yes 

thermal 

stress 

11A 

2kl 

lkjj 

1EX 

IE  n 

lax 

1  aii 

0.094 

Yes 

thermal 

stress 

11B 

2k  x 

lkjj 

1EX 

1E„ 

lax 

1«ii 

0.094 

Yes 

thermal 

stress 

12A 

lkx 

lkjj 

1EX 

lEu 

2a  j 

l°n 

0.094 

Yes 

thermal 

stress 

12B 

lkx 

lkrI 

1EX 

1E„ 

2*i 

la„ 

0.094 

Yes 

thermal 

stress 

•  Base  ■aterials  for  both  the  coating  layer  and  the  substrate  are 


Stellite  III. 


increase  the  principal  thermal  stress.  In  Figure  4.16,  the  material  of 
the  coating  layer  is  Stellite  III  for  all  cases.  Young's  modulus  of 
the  substrate  is  one-fifth  (case  10A),  one-half  (case  10B),  and  five 
times  (case  IOC)  that  of  Stellite  III.  This  figure  shows  that 
decreasing  Young's  modulus  in  the  substrate  will  result  in  increasing 
the  thermal  stress  field,  but  the  Influence  from  the  substrate  is 
weaker  than  from  the  coating  layer.  Figure  4.17  compares  the  effect  of 
Young's  modulus  on  the  thermal  stress  field  from  a  single  material  and 
from  a  layered  medium.  In  the  figure,  dashed  lines  represent  the  case 
of  a  single  material  with  a  cavity,  while  solid  lines  represent  the 
case  of  a  layered  medium  with  a  cavity.  From  this  figure,  we  observe 
that  thermal  stress  increases  linearly  in  proportion  to  Young's  modulus 
for  the  single  material  case.  For  the  layered  medium  case,  however, 
increasing  by  the  same  amount  Young's  modulus  in  the  coating  layer  will 
result  in  higher  thermal  stress  than  in  the  case  of  a  single  material. 
This  is  because  we  will  have  a  relative  softer  substrate  by  increasing 
Young's  modulus  in  the  coating  layer.  The  effects  from  thermal 
conductivity  and  the  coefficient  of  thermal  expansion  are  presented  in 
Figures  4.18  and  4.19.  These  effects  are  similar  to  those  found  in  the 
case  of  a  single  material  with  a  cavity. 

From  the  failure  specimen  for  the  case  of  a  single  material  with  no¬ 
cavity  (Figures  1.1  and  1.2),  we  observe  that  the  thermomechanical 
cracking  is  perpendicular  to  the  wear  surface.  However,  Ju  and  Liu 
C341  showed  that,  in  the  case  of  a  layered  medium  with  no-cavlty,  shear 
delamination  (cracking  is  parallel  to  the  wear  surface)  may  occur 
caused  by  the  changing  of  principal  directions  (larger  angle  of 
principal  direction),  therefore,  it  is  important  to  understand  what 
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4.17  Thermal  principal  stress  (cases  1A, SA.SC.7A, 9B> .  L=0.094,  and  n=0.06  above 


principal  stress  (cases  1A,12A,12B).  L=0.094,  and  n-0.06  above 


will  affect  principal  directions.  Figure  4.20  shows  the  effect  of  the 
ligament  thickness  (cavity  location)  on  principal  directions  at  the 
point  C=0.3  and  rpO.006  above  the  cavity  top  edge.  From  this  figure, 
we  observe  that  increasing  the  ligament  thickness  will  result  in  larger 
angle  of  principal  direction,  and  the  angle  of  principal  direction 
changes  drastically  when  L=0.094,  the  value  which  gives  the  maximun 
tensile  stress.  Figures  4.21  and  4.22  compare  the  effect  of  Young's 
modulus  on  principal  directions  for  the  case  of  a  single  material  with 
a  cavity  (dashed  line)  and  for  the  case  of  ^  layered  medium  with  a 
cavity  (solid  line).  These  two  figures  establish  that  decreasing 
Young's  modulus  in  the  coating  layer  (EJ>  or  increasing  Young's  modulus 
in  the  substrate  (E^)  will  increase  the  angle  of  principal  direction. 
Nevertheless,  changing  Young’s  modulus  in  the  case  of  a  single  material 
with  a  cavity  will  not  affect  the  principal  directions.  Figures  4.23 
and  4.24  illustrate  the  effects  of  the  thermal  conductivity  and  the 
coefficient  of  thermal  expansion  on  principal  directions  for  the  case 
of  a  single  material  with  a  cavity  (dashed  line)  and  the  case  of  a 
layered  medium  with  a  cavity  (solid  line).  From  these  two  figures,  one 
can  see  that  thermal  conductivity  and  the  coefficient  of  thermal 
expansion  will  not  influence  principal  directions  significantly. 
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direction  of  the  theraal  stress  field.  L-0 . 094 ,  and  n=0.06  above  the  top 


DISCUSSION  AMD  CONCLUSIONS 


The  Dresent  Investigation  demonstrates  the  effect  of  a  near¬ 
surface  rectangular  cavity  on  the  temperature  and  stress  fields  caused 
by  the  frictional  excitation  of  a  moving  asperity.  The  effects  of  a 
coating  layer  are  demonstrated  by  the  material  parameter  variations  in 
the  coating  layer  and  the  substrate,  including  changes  on  both  thermal 
and  mechanical  properties.  The  mathematical  model,  because  of  the 
geometry,  is  time  explicit.  Since  the  transient  solution  to  a  multiple- 
boundary  problem  is  always  complex,  numerical  solutions  become 
necessary  for  analyses  in  specific  cases.  In  the  present  problem,  it 
has  been  demonstrated  that:  <i>  the  transient  governing  differential 
equations  (2.1,  2.12,  and  2.26),  can  be  formulated  in  difference  forms; 

< 1 1 >  the  nonuniform  mesh  <t,q>,  which  must  be  employed  due  to  strong 
local  effects,  can  be  transformed  into  a  uniform  mesh  <£,n>;  <iii> 
boundary  conditions  in  temperature  and/or  heat  flux  can  be  expressed 
through  the  energy  balance  method,  thus  avoiding  the  singularity 
problem  at  the  cavity  corner;  (iv)  the  stress  singularity  at  each 
corner  of  the  cavity  can  be  taken  care  of  by  embedding  a  known  stress 
singularity  in  the  vicinity  of  the  corner;  (v)  the  numerical  solution 
can  be  tested  by  comparing  with  a  known  analytical  solution,  showing  a 
satisfactory  accuracy;  and  <vl)  the  numerical  scheme  can  be  extended  to 
compute  the  solution  for  other  geometries,  such  as  those  Including 
cracks  and  circular  cavities,  using  proper  coordinate  transformations. 

Like  most  numerical  solutions,  functional  relationships  can  not  be 


obtained  without  voluminous  computations.  However,  significant 
conclusions  can  be  reached  through  a  careful  selection  of  pertinent 
cases  for  the  numerical  results.  The  conclusions  for  the  present 
problem  are: 

Temperature  field 

1.  Because  of  the  discontinuity  in  heat  transfer  across  the  cavity, 
temperature  will  rise  higher  in  the  ligament  region  than  the  no-cavity 
case. 

2.  The  temperature  rise  is  inverse  to  the  ligament  volume, 
represented  by  the  ligament  thickness. 

3.  Increasing  the  thermal  conductivity  and  heat  capacity  of  the 
coating  layer,  will  decrease  the  surface  temperature. 

4.  When  the  coating/substrate  interface  is  at  the  ligament  depth, 
the  thermal  property  of  the  substrate  will  influence  the  temperature 
field  in  the  region  on  the  trailing  edge  of  the  asperity. 

5.  Because  of  the  necessary  heat  transfer  in  the  lateral  direction, 
the  heat  flux  will  be  at  a  large  oblique  angle  to  the  wear  surface.  In 
the  case  of  a  layered  medium  without  a  cavity,  the  near  surface  heat 
flux  at  the  critical  position  is  in  a  direction  approximately  90°  to 
the  wear  surface.  With  the  presence  of  a  cavity,  not  only  the 
magnitude  of  the  temperature  gradient  Increases,  but  also  the  direction 
of  the  temperature  gradient  is  rotated  to  a  more  oblique  angle  to  the 
wear  surface.  This  will  affect  principal  directions  in  the  thermal 
stress  field. 

Stress  field 

1.  In  the  governing  differential  equation,  the  small  order 
coefficients  of  the  dynamic  terms  would  have  required  an  extremely 


■sail  time  step  for  the  consideration  of  stability  and  truncation 
error.  This  difficulty  was  circumvented  by  using  the  perturbation 
method.  The  solutions  of  the  differential  equations  <4.12)  and  <4.13) 
of  the  various  perturbation  orders  are  well-behaved.  The 
magnification,  <ui+1/ui),  is  of  order  102.  Since  e<=M2>  is  of  the 
order  10”6,  each  perturbation  term  in  equations  <4.12)  and  <4.13)  is  of 
order  10~4  of  its  preceding  term.  Because  the  series  converges 
rapidly,  all  computations  are  deemed  adequate  by  using  only  one  term. 

2.  When  a  cavity  exists,  the  stress  state  that  determines  the 
failure  phenomenon  is  much  more  severe  than  in  the  no-cavity  case. 

This  will  lead  to  earlier  failure  of  the  mechanism. 

3.  The  mechanical  effect,  which  can  be  neglected  in  the  no-cavity 
case,  is  not  negligible  when  cavities  exist. 

4.  The  effects  of  the  mechanical  and  thermal  properties  on  the 
stress  field  are  consistent  with  those  obtained  in  the  no-cavity  case 
in  reference  C 32 1 .  These  effects  may  be  summarized  as  follows:  thermal 
stress  can  be  reduced  by  decreasing  Young's  modulus  in  the  coating 
layer,  increasing  Young's  modulus  in  the  substrate,  increasing  thermal 
conductivity  and  thermal  capacity  of  the  coating  layer,  and  decreasing 
the  coefficient  of  thermal  expansion  of  the  coating  layer. 

5.  For  a  thin  coated  medium,  the  cavity  location  and  the  material 
properties  matching  Especially  Young's  modulus)  will  influence  the 
principal  directions  of  the  thermal  stress  field.  When  the  angle  of 
principal  direction  becomes  larger,  shear  stress  at  the 
coating/substrate  Interface  becomes  dominant,  leading  toward 
delamination  of  the  coating. 

6.  The  location  of  the  cavity  influences  the  critical  depth  at 


which  the  thermal  tensile  stress  reaches  a  maximum.  When  the  cavity 
occurs  closer  to  the  wear  surface,  not  only  the  critical  depth  is 
reduced  but  also  a  higher  stress  results,  which  reaches  Its  maximum  at 
a  critical  ligament  thickness.  Further  reduction  of  the  ligament 
thickness  would  increase  the  ligament  temperature,  resulting  in  an 
extension  of  the  thermal  compressive  region  therein.  Correspondingly, 
the  thermal  tensile  stress  decreases  near  the  ligament  region.  The 
illustration  for  Stellite  III  shows  that  the  critical  thickness  is  at 
Lcr=0.094  for  both  cases  of  a  single  material  and  a  layered  medium  with 
a  thickness  of  approximately  40Z  of  the  critical  depth  of  the  no-cavity 
case.  For  the  normal  design  of  coating  thickness,  the  critical  depth 
of  the  specific  coating  material  can  function  as  a  guide.  However,  if 
cavities  are  either  unavoidable  or  too  expensive  to  control,  the  design 
thickness  should  avoid  the  critical  ligament  thickness. 
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APPENDIX  I 


nmCDDCTIGM  TO  THE  FINITE  DIFFERENCE  METHOD 


The  use  of  numerical  methods  for  solving  problems  is  a  result  of 
the  complexity  of  the  analytical  solutions  associated  with  practical 
engineering  problems.  Often  times,  analytical  solutions  are 
impossible.  In  engineering  problems,  factors  that  bring  about  the  use 
of  numerical  methods  are  complex  geometry,  nonlinearity,  nonuniform 
boundary  conditions,  time-dependent  boundary  conditions,  temperature- 
dependent  properties,  and  so  on.  In  some  cases,  analytical  solutions 
are  possible,  in  principle,  but  the  mechanics  of  obtaining  the  exact 
solution  may  be  much  more  difficult  than  the  task  of  solving  the 
problem  numerically.  For  example,  In  the  problem  of  finding  the  stress 
solution  of  a  composite  multilayered  body  with  nonhomogeneous  boundary 
conditions.  It  is  relatively  easy  to  set  up  the  differential  equations. 
The  solution,  however,  is  extremely  complex,  because  it  is  necessary  to 
deal  with  simultaneous  partial  differential  equations.  In  all  such 
cases  and  many  others,  if  one  is  equipped  with  the  knowledge  of 
numerical  methods  and  computer  programming,  the  required  solution  can 
be  successfully  obtained. 

Finite  difference  approximations  for  derivatives  were  already  in 
use  by  Euler  in  1768.  The  simplest  finite  difference  procedure  for 


dealing  with  the  problem  dx/dt=f(x>,  x(0)=a  is  obtained  by 
replacing  (dx/dt^-j  with  the  crude  approximation  <xn  -  xn«1>/At. 


* 


k 


k 


r  Vs 


that,  for  one-dimensional  systems,  the  finite  difference  approach  has 
been  deeply  Ingrained  in  computational  algorithms  for  quite  some  time. 
I-l  Finite  Difference  Approximation  of  Derivatives  Through 
Taylor 'a  Series 

The  derivative  of  a  function  at  a  given  point  can  be  represented  by 
a  finite  difference  approximation  using  a  Taylor  series  expansion  of 
the  function  about  that  point.  Let  f(x>  be  a  function  that  can  be 
expanded  in  a  Taylor  series.  Then  a  Taylor  series  expansion  of  the 
functions  f(x+h)  and  f(x-h)  about  x,  as  illustrated  in  Figure  1.1,  is 
given  by 


h2  h3 

f  (x+h)  =  f  ( x )  ♦  h  f*(x)  «•  —  f "  <x>  ♦  —  f*"(x>  ♦  ....  (I.l> 

h2  h3 

f(x-h)  =  f(x)  -  h  f '  (x)  +  —  ~  57  f*"<x>  ♦  ....  (1.2) 

where  primes  denote  derivatives  with  respect  to  x.  The  first-  and 
second-order  derivatives  f'<x>  and  f ' * ( x >  can  be  represented  in  the 
finite  difference  form  in  many  different  ways  by  utilizing  Taylor 
series  expansions  given  by  equations  (I.l>  and  (1.2)  as  now  described. 

First  Derivatives 

To  obtain  expressions  for  the  finite  difference  form  of  the  first- 
order  derivative  f*<x>,  equations  (1.1)  and  (1.2)  are  solved  for 
f'(x).  We,  respectively,  obtain 

f (x+h)  -  f ( x )  h  h2 

f  *  <  x )  =  - ^ - f”(x)  -  —  f'”(x)  ♦  (1.3) 


f(x>  -  f(x-h)  h  hz 

f '(X)  =  - r -  ♦  -r-  f  ”  (x>  -  —  f"’(x)  ♦  ...,  (1.4) 

n  4  b 
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Subtracting  equations  <I.1>  and  (1.2)  and  solving  for  f'(x)  we  obtain 


f(x+h)  -  f(x-h)  h2 

f(«l - 5J - r  f"«K>  *  .... 


(l.b) 


% 

9 

fi 


From  equations  (1.3)  to  (1.5),  the  following  approximations  can  be 
written  respectively  for  the  first  derivative  of  a  function  f(x)  about 
the  point  x. 


r  = 

i 


h 

f  - 

1  + 

Vi  + 

h 

f 

-  f 

_i±i- 

_ in 

0(h),  forward  difference 


2h 


(I.S) 


(1.7) 


(1.8) 


here  the  notation  0(h)  is  used  to  show  that  the  truncation  error 
involved  is  of  the  order  of  h;  similar’y  0(h2)  is  for  the  truncation 
error  of  the  order  of  h2,  and 


x  =  lh,  x+h  =  <i+l)h,  x-h  =  (i-l)h. 


etc, 


f ( x )  =  fit  f <x+h>  =  f1+j,  f(x-h>  =  fi_j,  etc. 


(1.9) 


( 1 . 10) 


V  V  ' 


as  illustrated  in  Figure  1.2.  We  note  that  forward  and  backward 
differences  are  accurate  to  the  order  h  whereas  the  central  difference 
expression  is  accurate  to  the  order  h2.  More  accurate  expressions  can 
be  obtained  for  the  forward  and  backward  difference  representation  of 
the  first-order  derivative  as  will  be  dlcussed  later . 
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Second  Derivatives 


We  now  proceed  to  the  finite  difference  representation  of  the 
second  derivative  f''(x)  of  a  function  f(x>  about  the  point  x.  To 
obtain  such  results  we  consider  a  Taylor  series  expansion  of  functions 
f(x+2h>  and  f(x-2h)  about  x  as 

4 

f (x+2h>  =  f ( x )  +  2h  f'(x)  +  2h2  f"(x)  +  —  h^  f"'(x>  +  ....  (1. 11) 

4 

f<x-2h>  =  f(x)  -  2h  f'(x>  ♦  2h2  f’’(x) - —  h3  f'”(x)  +  ....  (1.12) 

Eliainating  f'(x)  between  equations  (1.1)  and  (1. 11)  we  obtain 

f(x)  *  f(x+2h>  -  2f(x+h> 

f  "  <x)  =  - ^ - h  f"'(x>  +  ....  (1.13) 

Similarly,  eliminating  f'(x)  between  equations  (1.2)  and  (1.12),  we 
find 


f(x-2h)  +  f(x)  -  2f <x-h) 

f"(x)  =  - ^ -  +  h  +  •••• 


(1.14) 


Eliminating  f'(x)  between  equations  (I.l)  and  (1.2)  we  obtain 


f <x-h)  ♦  f (x+h)  -  2f(x)  1 

f '  '  (x)  =  - — -  -  —  h2  f""(x)  ♦  ....  (1.15) 


h2 


12 


Using  the  subscript  notation  defined  by  equations  (1.9)  and  (I. 10), 
various  forms  of  the  finite  difference  representation  of  the  second- 
order  derivative  f''(x)  about  the  point  x  given  by  equations  (1.13)  to 
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(I. IS)  are  written,  respectively,  as 


The  above  procedure  can  be  extended  to  obtain  more  accurate 
expressions  for  the  first  and  second  derivatives.  Such  expressions 
are  presented  in  reference  £47)  for  various  order  derivatives. 


The  derivative  of  a  function  in  non-uniform  spacing  (Figure  1.3)  can 
also  be  approximated  by  finite  difference  using  a  Taylor  series 
expansion.  A  summary  of  the  finite  difference  representation  of  the 
first-  and  second-order  derivatives  of  a  function  f(x)  in  non-uniform 
spacing  is  given  below 


f;  ■  -  e 


2ht  ♦  h2 


(h.+h~)  fi  + 


hi  +  h2 


hih2 


h2  (hi+h2 >  i+2  ’ 


i+i 

forward  difference 


(1.22) 


hi  <hi+h2> 


1-2 


hl  4  h2 
hjh2 


hx  +  2hz 


i-i 


backward  difference 


(1.23) 


f;  =  -h 


h2  -  hj 


i  (hi+h2>  i-i 


hih2 


fi  *  h2  (hi+h2)  fi+i  ’ 

central  difference 


(1.24) 


2  2  2 
^i  hi  (hi+h2>  ^l~i  hjhg  ^l  h2  (hi+h2>  ^i+i 

central  difference  '1.25) 

Using  equations  (1.22)  to  (I.2S)  is  very  cumbersome,  and  it  mav  lead 
to  loss  of  accuracy.  A  more  elegant  method,  general  coordinates 
transformation,  can  be  employed  in  the  non-uniform  mesh  to  avoid  these 
problems.  This  transformation  will  be  discussed  in  detail  later. 

1.2  Errors  Involved  in  lh»erlcal  Solutions 

In  numerical  solutions  using  the  method  of  finite  differences,  the 
partial  differential  equation  is  approximated  with  finite  difference 
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expressions  at  each  nodal  point,  and  as  a  result  the  solution  of  the 
differential  equation  is  transformed  to  the  solution  of  a  set  of 
algebraic  equations.  Wc  have  seen  that,  whenever  a  derivative  is 
approximated  by  finite  difference  using  a  Taylor  series  expansion,  an 
error  is  involved.  Such  an  error  is  called  the  truncation  error  or  the 
discretization  error.  These  errors  appear  because  a  continuous 
operator  such  as  the  first,  or  the  second-order  derivative,  is  replaced 
by  a  finite  difference  approximation.  In  addition,  numerical 
calculations  are  carried  out  only  to  a  finite  number  of  decimal  places 
or  significant  figures;  as  a  result,  at  each  step  in  the  calculation, 
some  error  is  introduced  due  to  this  rounding-off,  called  the  round-off 
error. 

Clearly,  if  the  finite  difference  approximation  is  made  by  using 
formulas  having  truncation  errors  of  high  order,  the  truncation  error 
at  each  step  is  minimum.  Also,  by  decreasing  the  step  size,  the 
truncation  error  is  reduced  for  each  step;  however,  a  limit  also  is 
reached  at  which  further  reduction  in  step  size  increases  the  total 
number  of  calculations  and  as  a  result  the  round-off  error  may  become 
dominant. 

Ideally,  if  it  were  possible  to  carry  out  the  finite  difference 
calculations  with  extremely  small  steps  and  to  perform  the  calculations 
to  an  infinite  number  of  decimal  places,  the  resulting  solution  would 
be  exact.  However,  due  to  the  cumulative  effects  of  the  rounding  off 
error  and  the  discretization  errors,  the  solutions  obtained  by  the 
finite  difference  method  is  expected  to  deviate  from  the  exact  result; 
therefore,  the  solution  computed  is  the  numerical  solution  but  not  the 
exact  result.  It  is  very  difficult  to  determine  the  cumulative 


departure  of  the  numerical  solution  from  the  exact  result  due  to  the 
cumulative  effects  of  such  errors.  Comparison  of  numerical  solutions 
with  exact  analytic  solutions  reveals  that,  for  most  cases,  the  results 
are  very  close  Indeed.  After  some  experience  with  different  methods 
and  different  step  sizes,  a  suitable  combination  can  be  chosen  for  the 
numerical  solution  of  a  given  problem. 

1.3  Time  Dependent  Problem 

The  stability  consideration  plays  an  important  role  in  the  finite 
difference  solution  of  time  dependent  problems.  There  are  several 
schemes  available  to  express  the  time  dependent  problems  in  finite 
difference  form.  Each  of  these  differencing  schemes  has  its 
advantages  and  limitations.  We  now  discuss  some  of  them  by  using  the 
one-dimensional  time-dependent  heat  conduction  equation  as  examples. 
Explicit  Method 

The  one-dimensional,  time-dependent  heat  conduction  equation  for  a 
finite  region  0  <  x  <  L  is 


3T  32T 
3t  =  k  3x2  * 


(1.26) 


32T  dT 

If  ar,d  3^  are  replaced  by  the  central  and  forward  differences, 

respectively,  and  using  a  uniform  mesh  size  Ax  in  the  x  domain  and  At 
for  the  time  step,  equation  (1.26)  can  be  rewritten  in  the  finite- 
difference  form  as 


T(i,n+1>  -  T<i,n>  T<i-l,n>  -  2T(l,n>  +  T(i+l,n) 
_  =  k  Ax 2 


(1.27) 


with  a  truncation  error  of  0(At >+0(Ax2 ) .  Where  T(x,t>  is 


represented  by  T(x,t>  =  T(iAx,nAt)  =  T(i,n>.  Solving  equation  <1. 27) 
for  T(i,n+1)  one  obtains 

T<i,n+1>  =  rT<i-l,n)  <l-2r  )T(  i  ,r.)  +  rT(i+l,n>,  <1.28) 

kAt 

where  r  =  • 

The  equation  is  called  the  explicit  form  because  the  unknown 
temperature  T<i,n+1)  at  the  time  step  (n+i)  can  be  directly  determined 
from  the  temperatures  T<i-l,n>,  T(i,n>,  and  T<i+l,n>  at  the  previous 
time  step.  The  explicit  scheme  provides  a  relatively  straightforward 
expression  for  the  determination  of  the  unknown  T(i,n+1>.  The  only 
disadvantage  of  this  method  is  that  once  k  and  Ax  are  fixed,  there  is  a 
maximum  permissible  step  size  At,  which,  by  instability  considerations 
should  not  be  exceeded.  For  example,  when  the  boundary  conditions  at 
x=0  and  x=L  are  both  of  the  first  kind  (i.e.,  specified  temperature), 
the  restriction  imposed  on  the  parameter  r  is 

1 

0  <  r  <  -  .  (1.29) 


That  is,  for  given  values  of  k  and  Ax,  if  the  time  step  At  exceeds  the 
limit  Imposed  by  the  above  criteria,  the  numerical  calculations  become 
unstable,  as  a  result  of  the  amplification  of  errors.  Figure  1.4 
illustrates  what  happens  to  the  numerical  calculations  when  the  above 
stability  criteria  is  violated.  In  this  figure,  the  numerical 


calculations  performed  with  a  time  step  satisfying  the  condition 
S  1 

r  =  ^  is  in  good  agreement  with  the  exact  solution;  whereas  the 


numerical  solution  of  the  same  problem  with  slightly  larger  time  step, 
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Finite  difference  solution 


1.4  Effects  of  parameter  r»*At/<Ax>2  on  the 
stability  of  finite  difference  solution. 


'Z/V/V/V 


which  violates  the  above  stability  criteria,  i.e.,  r=->-  >  results 


an  unstable  solution. 


Implicit  He t hod 


The  explicit  method  discussed  above  is  simple  computationally,  but 

very  small  time  step  should  be  used  because  of  stability 

considerations.  Therefore,  a  prohibitively  large  number  of  time  steps 

may  be  required  if  solutions  are  to  be  computed  over  a  large  period  of 

time.  It  is  for  this  reason  that  other  finite  difference  forms,  found 

to  be  less  restrictive  to  the  size  of  time  step  At,  have  been 

developed.  One  such  scheme  is  the  fully  Implicit  method.  Ue 

illustrate  this  method  by  considering  the  finite  difference 

representation  of  the  heat  conduction  equation  (1.26).  The  partial 
32T 

derivative  ^5  ls  represented  in  finite  difference  form  using  the 

3T 

central  difference  formula,  whereas  the  time  derivative  —  is 

dt 

represented  in  the  finite  difference  form  using  the  backward  difference 
expression. 

Then,  the  finite  difference  form  of  equation  (I.26>  becomes 


T<i,n+1>  -  T(i,n>  T(i-l,n+l)  -  2T<i,n«-l>  ♦  T(i*l,n+1) 

- - -  - — -  ,  ,1.30, 


This  is  called  an  implicit  form  of  the  finite  difference  representation, 
because  to  determine  the  unknowns  T(i,n+1),  a  set  of  simultaneous 
algebraic  equations  are  to  be  solved.  The  advantage  of  the  implicit 
method  is  that  it  is  stable  for  all  sizes  of  time  step  At.  Thus,  there 
is  no  size  restriction  on  At.  The  only  size  restriction  on  At  is  due 
to  the  consideration  of  the  truncation  error. 

The  truncation  errors  for  both  explicit  and  implicit  forms  of  the 
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finite  difference  representations  of  the  heat  conduction  equation  Is 
of  the  order  (Ax2)+(At).  But  the  actual  accumulated  error  in  both 


methods  need  not  be  the  same.  Depending  on  the  nature  of  the  problem, 


one  of  the  methods  may  be  preferred  to  the  other. 


Crank-Nicolson  Method 


Crank  and  Nicolson  [471  suggested  a  modified  implicit  method.  In 


this  method,  the  heat  conduction  equation  (1.26)  is  represented  in 


finite  difference  form  by  taking  the  arithmetic  average  of  the  right- 


hand  sides  of  the  explicit  form  (1.26)  and  the  implicit  form  (1.30). 


Then,  equation  (1.26)  becomes 


T< i , n+1 )  -  T< i ,n) 


k  T(i-l,n+l>  -  2T< i , n+1 >  +  T(i,n+1> 

2  C  Ax2  + 


T(i-l.n)  -  2T(i,n)  +  T(i+l.n) 

+  - Z? -  ]- 


(1.31) 


The  advantage  of  this  mehtod  is  that,  for  given  values  of  the  space  and 


time  steps  Ax  and  At,  the  resulting  solution  involves  less  truncation 


error  due  to  At  than  the  explicit  and  the  implicit  forms  discussed 


above.  On  the  other  hand  the  Crank-Nicolson  form  involves  additional 


computation. 


To  provide  a  better  insight  to  the  physical  significance  of  the 


Crank-Nicolson  representation,  equation  (1.31)  can  be  written  in  a  more 


general  form  by  taking  a  weighted  average  of  the  two  terms  in  the 


brackets 


T(i,n+1)  -  T(i,n>  T<i-i.n+l)  -  2T(i,n+l)  ♦  T(i+l,n+l> 

- - -  -  k  l  r,  - ^ -  ♦ 


where  0  <  n  <  1  is  called  the  degree  of  implicitness.  Clearly, 

equation  (1.32)  reduces  to  the  explicit  form  given  by  equation  (1.27) 

for  rj=0,  to  the  implicit  form  given  by  equation  (1.30)  for  rpl  and  the 

1 

Crank-Nicolson  form  (1.31)  for  q=-  . 

Alternating-Direction  Implicit  Method 

The  Implicit  methods  discussed  above  are  advantageous  to  us  because 
of  the  superior  stability  properties.  On  the  other  hand,  becuase  a 
large  number  of  simultaneous  equations  need  to  be  solved  at  each  time 
step,  the  computational  problems  become  enormous  when  they  are  applied 
to  the  solution  of  time  dependent  problems  involving  two  or  three  space 
dimensions.  For  example,  for  a  three-dimensional  problem  with  N 
interior  nodal  points  in  each  direction,  there  are  a  total  of  N3  nodes, 
hence  N3  x  N3  matrix  equations  must  be  solved  for  each  time  increment. 

The  alternating-direction  implicit  (A.D.I.)  method  introduced  by 
Peaceman  and  Rachford  [481,  provides  an  efficient  method  for  solving 
problems  involving  large  number  of  nodes.  To  illustrate  the  procedure, 
a  two-dimensional,  time  dependent  heat  conduction  equation  is 
considered 

32T  32T  1  8T 

*  dp  =  K  at  ’  in  the  region  R  (1.33) 

To  represent  the  space  derivatives  in  finite  difference  form,  the  cent¬ 
ral  difference  formula  is  used  with  an  implicit  and  explicit  difference 

32T  32T  32T 

approximation  alternatively  for  ®n<*  •  For  example,  if 


32T 

is  represented  in  the  implicit  form,  the  derivative  g^g  is 
represented  by  an  explicit  approximation.  Then,  the  finite  difference 
form  of  equation  (1.33),  to  proceed  the  solution  from  the  <n>th  step 
to  (n+l)th  step,  becomes 


1  T(i,j, n»l>  ~  T(i,j,n>  T(i-l,j,n+D  -  2T(i,j,n+l>  +  T(  i+1 ,  J ,  n+1 ) 
k  At  Ax2 


T(i,j-l,n>  -  2T ( i , j , n )  +  T(i,j+l,n) 
Ay2 


(1.34) 


The  finite  difference  form  of  equation  (1.33),  to  proceed  the 

solution  from  the  (n+l)th  step  to  the  (n+2)th  step,  is  written  using 

32T  32T 

an  explicit  form  for  g^j  and  implicit  form  for  g^  as 

1  T<  1 , j , n+2)-T( i , j  ,n+l )  T( i-1 , J ,n+l )  -  2T(l,j,n+l>  +  T(i+1, j,n+l> 
k  At  "  Ax2 

T(i, 1-1, n+2)  -  2T(i,1,n+2)  ♦  T(i,J+l,n+2) 

+  - - — 5 -  .  (1.35) 

Ay2 

The  procedure  is  repeated  alternately  in  the  subsequent  time  steps. 

The  advantage  of  the  A.D.I.  method  over  the  implicit  method  results 
from  the  fact  that  it  reduces  the  number  of  equations  to  be  solved 
simultaneously  for  each  time  step.  Consider  for  example  a  two- 
dimensional,  time  dependent  problem  with  N  internal  nodes  along  the  x- 
axis  and  N  nodes  along  the  y-axis.  The  A.D.I.  method  requires  the 
solution  of  N  simultaneous  equation  N  times  for  each  time  step,  whereas 
the  implicit  method  requires  the  solution  of  N2  equation  at  each  time 
step.  There  are  other  methods,  for  example,  the  alternating  direction 
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,1 

i 


explicit  (A.D.E)  method  [49],  the  Douglas-Rachford  implicit  scheme 


C481,...,  etc.  The  reader  should  consult  these  references  for  further 


discussion  of  these  methods. 


APPENDIX  II 


WI FORM  MESH  AND  GENERAL  COORDINATES  TRANSFORMATION 


The  solution  of  a  system  of  partial  differential  equations  can  be 


greatly  simplified  by  a  well-constructed  grid.  On  the  other  hand,  a 


grid  which  is  not  well  suited  to  the  problem  can  lead  to  an  unsatis¬ 


factory  result.  In  some  applications,  improper  choice  of  grid  point 


locations  can  lead  to  an  apparent  instability  or  lack  of  convergence. 


For  many  applications,  a  non-uniform  mesh  must  be  used  in  order  to 


obtain  an  accurate  solution  and  to  save  computing  time.  One  can  solve 


the  problem  in  the  physical  plane  (original  plane)  by  applying  the 


difference  formulas  on  the  non-uniform  mesh  directly,  or  transform  the 


non-uniform  mesh  to  a  uniform  mesh  and  solve  the  problem  in  the 


computational  plane.  Generally,  the  coordinate  transformation  gives  a 


more  accurate  solution  than  mesh  changes. 


II. 1  Non-Uniform  Mesh 


The  simplest  variation  of  the  rectangular  mesh  system  is  obtained 


by  simply  changing  the  mesh  spacing  in  one  direction  at  some  point. 


This  would  be  done  for  the  purpose  of  obtaining  higher  resolution  (and 


hopefully  higher  accuracy)  in  some  region  where  the  gradients  were 


expected  to  change  rapidly.  To  illustrate  this  technique,  we  consider 


the  obvious  method  of  changing  from  Axs  to  Ax2  between  node  points  at 


some  node  i=m,  as  shown  in  Figure  II.l. 


Expanding  a  function  in  a  Taylor  series  forward  and  backward  from  i=m 


gives 


V 


f  .  =  f  ♦  r-  Ax2  +  -  Ax2  +  Z  Ax23  ♦  0(Ax2«),  (II.l) 

»+ 1  m  ox  m  2  2  3x2  m  2  6  oxJ  m  2  2 


3f  1  32f  1  33f 

f  =  f  -  r -  Ax.+in  Ax  2  -  -  —  Ax  3  +  CXAx.**),  <11. 2) 

m~i  a  Ox  m  i  2  3x2  m  l  6  3x3  m  1  1 


The  expression  for  —  is  obtained  by  subtracting  equation  (I I. 2) 
dx  m 


from  < I I . 1 > 


3f  1  32f 

f  -  f  =  r-  (Ax,  +  Ax,)  +  -  r-s  (Ax  2  -  Ax ,2)  +  0(Ax3) , ( II . 3 ) 
m+i  m~i  3x  m  1  2  2  3x2  m  2  l 


where  by  0(Ax3)  we  mean  the  largest  of  (XAXj3)  or  0<Ax23).  Solving  for 
3f  I 

8lves 


»f I  ~  i  »2n  Ax  ,2  -  Ax ,2 


3x  ’  Axt  ♦  Ax2  2  3x2[m  Axt  +  Ax2 


+  0(Ax2), 


(II. 4) 


This  means  that  the  form 


f  ~  f 

■♦l _ nr 


dx | m  Axj  ♦  Ax 2  ’ 

is  second-order  accurate  only  if 


(II. S) 


Ax, 2  _  Ax ,2 

o  i  a — rr  3  <  o(Ax  2). 

Axj  ♦  Ax 2  -  * 


(11.6) 


Note  that,  for  Ax2  very  small,  the  accuracy  at  m  deteriorates  to  first 
c^der  in  Ax,. 

The  expression  for  the  second  derivative  is  obtained  by  multiply¬ 
ing  equation  (II. 2>  by  s2=< Ax,/Ax j >2  and  adding  the  result  to  (II. 1). 


f  .  ♦  <1+S2)f  +  S2f  =  — |  Ax0(l~S)  +  T~r  Ax.,2  + 

OX^  m  * 


m  m 


-1  3x  m  2 


i  a3f 

+  Z  a~a  Ax.,2  (Axy  -  Ax.  >  ♦  CMAx*1), 
b  oXJ  m  ^  ^  1 


(II. 7) 


02f  f  +  <l+s2)f  +  s2f  1  -  s 

°  ‘  _  m+  l _ m _ m~i  _  __  f  * _ ®  ■> 

3x2  m  Ax„2  3x  m  1  Ax,  ) 


+  0[<Ax2  -  Ax t ) ,  Ax2I 


(II. 8) 


The  resulting  expression  now  requires  s=0(l-Ax,2)  just  to  be 


first-order  accurate  at  i=m. 


It  is  clear  from  the  above  equations  that,  unless  the  mesh  spacing 


is  changed  slowly,  the  formal  truncation  error  is  actually 


degraded,  rather  than  improved. 


I 1.2  Coordinates  Transformations 


Early  work  using  finite  difference  methods  was  restricted  to 


problems  where  suitable  coordinate  systems  could  be  selected  in  order 


to  solve  the  governing  equations  in  that  base  system.  As  experience 


in  computing  complex  problem  was  gained,  general  mappings  were 


employed  to  transform  the  physical  plane  into  a  computational  domain. 


Numerous  advantages  accrue  when  this  procedure  is  followed.  For 


example,  when  the  untransformed  equations  are  differenced  in  the 


expanding  mesh,  the  result  is  a  deterioration  of  formal  accuracy,  as 


we  have  seen;  but  the  transformed  equations  may  be  differenced  in  a 


regular  mesh  (such  as  constant  Ax,  Ay)  with  no  deterioration  in  the 


formal  order  of  truncation  error,  except  that  it  will  now  be  0(Ay2) 


rather  than  0(Ay2),  also,  the  body  surface  can  be  selected  as  a 


boundary  in  the  computational  plane  permitting  easy  application  of 
surface  boundary  conditions.  In  general,  transformations  are  used 
which  lead  to  a  uniform  space  grid  in  the  computational  plane  while 
points  in  physical  space  may  be  unequally  spaced.  This  situation  is 
shown  in  Figure  11.2.  When  this  procedure  is  used,  it  is  necessary  to 
include  the  derivatives  of  the  mapping  in  the  differential  equation. 
II. 2.1  Simple  Transformations 

In  this  section,  simple  independent  variable  transformations  are 
used  to  illustrate  how  the  governing  equations  are  transformed.  As  a 
first  example,  the  problem  of  clustering  grids  near  a  wall  is 
considered.  Figure  I I. 3a  shows  a  mesh  above  a  flat  plate  in  which 
grid  points  are  clustered  near  the  plate  in  the  normal  direction  (y). 
While  the  spacing  is  not  uniform  in  the  y  direction,  it  is  convenient 
to  apply  a  transformation  to  the  y  coordinate  so  that  the  governing 
equations  can  be  solved  on  a  uniformly  spaced  grid  in  the 
computational  plane  <x,y)  as  seen  in  Figure  II. 3b.  A  suitable 
transformation  for  a  two-dimensional  problem  is  given  by 
Transformation  Is 


(II. 9) 


y  =  i  - 


ln{C3  +  1  -  <y/h) 1/Cg  -  1  ♦  (y/h)]} 
InCfg  ♦  l>/(3  -  1>] 


1  <  3  <  »  (II. 10) 


This  stretching  transformation  clusters  more  points  near  y=0  as  the 
stretching  parameter  3  approaches  1. 

In  order  to  apply  this  transformation  to  the  governing  equations, 
the  following  partial  derivatives  are  formed 


-  -  _ 


Axt  I  4*2  I  ^2 


Fifire  II. 1  ■omnlfora  spacing. 


<b)  computational  plane 


Flgwe  II. 2  Happing  to  cosput 


ft 


3x  "  3x  3x  dx  3y  ’ 


(II. 11 ) 


3_  3x  3_  3y  3_ 
3y  ”  3y  3x  3y  3y  ’ 


(11.12) 


where 


dx  3y  3x 

3x  =  1  ’  3x  =  °  *  3y  =  °  ’ 


3y _ 2p _ 

and  3y  '  h(3z  -  [1  -  (y/h)]2}lnC<3  +  l)/(3  -  1)1 


As  a  result,  the  partial  derivatives  simplify  to 


3  3 

3x  =  35  ’ 


(11.13) 


3_  3y  3_ 
3y  "  *3y )  3y  ’ 


(11.14) 


If  we  now  apply  this  transformation  to  the  following  equation 


3u  3u 
3x  +  3y  =  °’ 


(II. IS) 


the  following  transformed  equation  is  obtained 


3u  3y  3u 

—  ♦  <  — >  rz  =  0. 
3x  3y  3y 


(11.16) 


This  transformed  equation  can  now  be  differenced  on  the  uniformly 


spaced  grid  in  the  computational  plane.  We  note  that  the  expression 


for  the  derivative  —  contains  y  so  that  we  must  be  able  to  express  y 
as  a  function  of  y.  This  is  refered  to  as  the  inverse  of  the 
transformation.  For  the  present  transformation,  given  by  equations 
<11. 9)  and  <11. 10),  the  inverse  can  be  readily  found  as 


x  =  x  , 


<3  ♦  i>  -  <3  -  i>{c<3  +  i>/<3  - 
c <3  +  i>/<3  -  i)Ti*v  +  i 


(11.18) 


Transformation  2: 


x  =  x  , 


<11.19) 


y  =  a  +  <1  - 


ln((($+[y(2<x+H/h]-2a]/{$-[y( 2a+l  )/hJ*Zaj) 
a>  ln[ <3+1 )/<3~l > 1  •  <1I,20> 


For  this  transformation,  if  a-0  the  mesh  will  be  refined  near  y=h 

1 

only,  whereas,  if  o=~  the  mesh  will  be  refined  equally  near  y=0 
(Figure  II. 4).  It  has  been  shown  that  the  stretching  parameter  3  is 
related  (approximately)  to  the  nondlmensional  thickness  <6/h>  by 


<  1  -  6/ h)1' i  ’ 


0  <  6/h  <  1 


<11. 21) 


where  h  is  the  height  of  the  mesh.  For  the  transformation  given  by 

3y 

equations  (11.19)  and  (11.20),  the  derivative  —  is 


3y  23  < 1  -  o ) <  2a  +  1 > 

3y  =  h{32  -  Cy< 2a  +  l)/h  -  2a]2}lnt<3  ♦  l>/<3  -  1)]  ’ 


(11.22) 


and  the  inverse  transfornat ion  becomes 


x  =  x  , 


(11.23) 


<0  +  2a ) C  ( 0  ♦  1 )/ <0  -  in<y-a)/(l-a)  -  g  +  2a 
y  =  h  72a  +  1)11  +~770  ♦  1 )/ ( 0  -  i ) ] <y-a)/( 1-a) ) 


(11.24) 


A  useful  transformation  for  refining  the  mesh  above  some  interior 
point  yc  (see  Figure  II. 5)  is  given  by 
Transformation  3: 


x  =  x  , 


(11.25) 


1  y 

y  =  B  +  -  sinh-1  t ( — 

1  y  c 


-  l)sinh(tB>]  , 


(11.26) 


where 


1  ,  1  ♦  (eT  -  1  Myc/h> 

B  =  2t  ln^  1  *  (e-T  -  1 ) (y_/h>  >  ’ 


0  <  T  <  ® 


(11.27) 


In  this  transformation,  x  is  the  stretching  parameter  which  varies 

from  zero  (no  stretching)  to  large  values  which  produce  the  most 

By 

refinement  near  y=y  .  The  metric  —  and  y  become 


By  slnh(TB) 

By  ”  TyTTl- T- rTyTyT^^l^sTnh^TrBTTTTi  * 


(11.28) 


y  =  yc  ♦ 


sinh[T(y  -  B>3 
sinh(xB) 


(11.29) 


For  our  final  transformation,  a  simple  transformation  which  can  be 


<»)  physical  plane 


7 


(b)  computational  plane 

Figure  II. S  Grid  clustering  near  an  interior  point. 


I I. 2. 2  Generalized  Transformation 

In  the  preceding  section,  the  sisple  independent  variable 
transformations  which  make  it  possible  to  solve  the  governing  equations 
on  a  uniformly  spaced  computational  grid  were  examined.  Let  us  now 
consider  a  completely  general  transformation  of  the  form 


t  =  t<x,y,z> , 


<11.35) 


n  =  n<x,y,z) , 


(11.36) 


C  =  t(x,y,z) , 


(11.37) 


which  can  be  used  to  transform  the  governing  equations  from  the 


physical  domain  (x,y,z>  to  the  computational  domain  (t,q,t>.  Using 


the  chain  rule  of  partial  differentiation,  the  partial  derivatives 


become 


3x  t'x  at  +  n>t  an  +  Cx  at  ’ 


(11.38) 


3y  3t  *  nV  3q  *  '’V  3(;  ’ 


(11.39) 


3z  "  3t  +  nz  an  +  Cz  at  * 


(11.40) 


The  derivatives  <tx,  nx.  tx,  ny.  ty,  tz,  nz.  tz>  appearing  in  these 


equations  can  be  determined  in  the  following  manner.  We  first  write 


the  differential  expressions 


dt  =  Cxdx  +  tydy  ♦  tzdz  , 


(11.41) 


dq  =  qxd*  ♦  nydy  +  nzdz  » 


(11.42) 


dt  =  txdx  +  tydy  ♦  tzdz  . 


(11.43) 


I  <1  «  a’* 


which  in  matrix  form  become 


Sy  Sz 


In  a  like  manner  we  can  write 


Therefore, 


/*x  ty 

I  nx  ny  n2 

Ux  cy  c2 


xc  xn  *c  \  -1 
yn  yc 

\  zn  \  I 


(11.44) 


(11.45) 


Vc  "  yczn 


-(x  z^  -  x  z  ) 

n  c  c  n 


"<Vc  -  yQzc>  x^  -  x^zj. 


Vn  ‘  Vc 


-(x_z  -  X  zr  > 

C  n  n  t 


Vc  *  Vn 


*(Vc  •  xcV 


(11.46) 


Vn  '  Vt 


Thus,  the  derivatives  are: 


=  J<y_z_  "  y_z  >, 


(11.47) 


Ey  ‘  'J<xnzc  ‘  xcV' 


e*  *  J<xnyc  '  xcyn’' 


nx  =  '  vtzc>. 


Hy  =  J(X{_2{.  -  X^Zj.  )  , 


n2  =  -J(x  y^  -  Xj-Yj-*. 


C„  =  J<yr.z  -  y  z_ ) , 

x  7s  n  n  t 


cv  =  -J  <  xrz  -  x  z_  > , 

y  t  n  n  S 


C-  =  J(x_.y  -  x  yr  > . 

2  Z  n  n  ( 


where  J  is  the  Jacobian  of  the  transformation 


a<s,n.o 

J  =  3(x,y,z)  =  nx  nv  n* 


*x 


^x  Cy  Cz 


which  can  be  evaluated  in  the  following  manner 


l\  \ \\  M 

.  8<x,y,z)  .  I  \ 

j  =  i/j-i  =  i/(  ^^7 )  =  I  yF  >v  I 


K  7n  7c 


z_  z 
t:  n  c 


i^^&ssise 


VAINV 


(11.48) 


(11.49) 


(II. SO) 


(II. SI) 


(11.52) 


(I I. S3) 


(II. 54) 


(11.55) 


(11.56) 


(11.57) 


<11. 58) 


=  1/txC<ynWn)  " 


xn<ySzc-yc2V 


+  x  <yrz  -y  zr  )  ] . 

c  C  n  n  C 


The  coordinate  derivatives  can  be  readily  determined  if  analytical 
expressions  are  available  for  the  inverse  of  the  transformation: 


x  =  x<C,n.O, 


(11.59) 


y  =  y<S,n.c>, 


(I I. SO) 


z  =  z<c,n.c>. 


<11.61 ) 


For  a  two-dimensional  problem,  the  derivatives  of  a  function 
f(x,y)  in  the  transformed  plane  (£,n>  are: 


f  =  <y  f  -  y  f  )/J  , 
x  n  t  Cl 


(11.62) 


fv  =  <x_f  -  x  f_)/J  , 
y  £  n  i  t 


(11.63) 


’  ‘Wee 


‘  ‘Vn  *  WyE-)  *  W-in^Vn  '  Ve’"3'1"'65’ 


f__,  =  (x  2frr  ~  2x  x  f  +  x  2f  )/J2  + 

‘yy  n  ££  £  n  £n  £  nn 


+  [<x  2yrr  -  2x_x  yr  +  x  2y  ><x  f  -  x  f  )  ♦ 

n  ££  £  n  £n  £  nn  n  £  £  n 


+  (x  2xrr  -  2x  x  x  +  x  2x  Hyrf  -  y  fr)]/J3.  (11.66) 

n  «  £  n  £n  £  nn  £  n  n  £ 


The  Laplacian  is  given  by 


V2f  =  (AifFr~2A2frn+A3fnn>/J2  *  C  < ' 


'  'VrfVE1  *  IA>y«W^A3V‘\VxEV1/J3’  <U  t1' 


or 


V2f  = 


(Ajf££  -  2A2f^  + 


A„f  +  A«f 

3  nn  q  r 


+  A5f^)/J2. 


(11.68) 


where 


(11.69) 


x,_x 

£  n 


+  y£yn 


(11.70) 


(11.71) 


A,. 


< Y„A,  -  x^A,)/J  , 


(11.72) 


A.-  =  (x  A,  -  y  A,  )/J 

5  n  7  n  6 


(11.73) 


A,;  =  A ,xrr  -  2A,xr  +  A~x 

6  1  2  £n  3  nn 


(11.74) 


A7  =  -  2A2y(;n  +  A3Vr 


(11.7b) 


Likewise,  the  Gradient  is  given  by 


Vf  =  i(y  fr  -  yrf  )e,  ♦  (x_f  -  x  f_)e,l/J 

nL  t  I]  1  f  r>  n  F  2 


Li  It 


(11.76) 


and  the  Divergence  is  expressed  by 


V.F  =  ty^F^  -  y(;<F1>n  *  x^F^  -  x^F^l/J 


(11.77) 


where  F  =  Ftet  ♦  F2e2  . 

Finally,  the  Curl  may  be  written  as 


(11.78) 


Curl  F  =  e,.Cy  (F->r  -  yr(F2>  -  xr(F.>  ♦  x  (F.J-.l/J  ,  (11.79) 

J  n  *  t  t  1  n  n  l 


where  J  =  x„y  -  x  y_  is  the  Jacobian  of  the  transformation,  and  the 

C  n  n  C 

subscripts  (x,y,C,n>  denotes  partial  derivatives  in  those  coordinates, 
respectively. 

For  cases  where  the  transformation  is  the  direct  result  of  a  grid 
generation  scheme,  the  metrics  can  be  computed  numerically  using 
central  difference  in  the  computational  plane. 

The  general  coordinate  transformation  can  be  employed  to 
transform  very  complicated  curvilinear  coordinates  to  simple 
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o 


rectangular  coordinates.  Some  examples  which  transform  the  problem  in 
physical  domain  to  computational  domain  are  given  in  Figures  11.7,8,9. 
The  detail  and  more  complicated  transformations  are  found  in  [49, SO]. 
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APPENDIX  III 


ENERGY  BALANCE  METHOD 


The  subjects  of  thermodynamics  and  heat  transfer  are  highly 
complementary.  For  example,  heat  transfer  is  an  extension  of 
thermodynamics  in  that  it  considers  the  rate  at  which  energy  is 
transported.  Moreover,  in  many  heat  transfer  analyses  the  first  law  of 
thermodynamics  (the  law  of  conservation  of  energy)  plays  an  important 
role.  In  our  application  of  the  conservation  laws,  we  first  need  to 
identify  the  control  volume,  a  fixed  region  of  space  bounded  by  a 
control  surface  through  which  energy  and  matter  may  pass.  With  respect 
to  a  control  volume,  a  form  of  the  energy  conservation  requirement  that 
is  most  useful  for  heat  transfer  analyses  may  be  stated  as  follows. 

"The  rate  at  which  thermal  and  mechanical  energy  enters  a  control 
volume  minus  the  rate  at  which  this  energy  leaves  the  control  volume 
must  equal  the  rate  at  which  this  energy  is  stored  in  the  control 
volume . " 

If  the  inflow  of  energy  exceeds  the  outflow,  there  will  be  an 
increase  in  the  amount  of  energy  stored  (or  accumulated)  in  the  control 
volume;  if  the  converse  is  true,  there  will  be  a  decrease  in  energy 
storage.  If  the  inflow  of  energy  equals  the  outflow,  then  a  steady 
state  condition  must  prevail  in  which  there  will  be  no  change  in  the 
amount  of  energy  stored  in  the  control  volume.  For  the  heat  conduction 
problem,  the  inflow  and  outflow  energy  of  the  control  volume  is  the 
heat  flux  Q=-kVT.  The  law  of  conservation  of  energy  is  the  key  concept 
of  the  energy  balance  method. 
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III.l  Energy  Balance  on  the  Cavity  Boundaries  (Figure  3.3) 

On  face  AB  (Figure  III.l).  Let  us  examine  an  arbitrary  node 
P(i,j),  on  the  face  AB,  surrounded  by  nodes  R,  S,  and  W.  When  heat 
conduction  from  the  boundary  nodes  R  or  W  to  P  is  examined,  we  observe 
that  the  area  available  for  heat  flow  is  only  I (Ay/2) -l],  although  the 
distance  across  which  heat  is  conducted  is  still  Ax.  Thus  we  have 

Tr(i,j-1)  -  Tr(i,j> 

Qs-»p  -  Rj(Ax)  ^  ,  (III.l) 

Tr(i-l,j>  -  T1 ( i , j  > 

Qw-»p  ~  ^  ^ ( Ay / 2 )  ^  ,  (III. 2) 


0R->P  =  k^Ay/2) 


TI(i-»l, j)  -  TI(i,j) 
Ax 


(III. 3) 


The  rate  of  change  of  internal  energy  Up  in  the  time  interval  At  at 
P(i,j)  is 


UP=  PjC^Ax-Ay/2)  - — 


TI(i,j,n)  -  TI(i,j,n-l) 


n=l  (III. 4a) 


UP=  pjci(AxAy/2) 


3Tr(i, J,n)-4TI(i, J,n-l)+TI(i, j,n-2) 
2At 


n>l  (III. 4b) 


From  conservation  of  energy,  Up  =  Qsum,  we  obtain 


T1 < i , j , n>  =  T1 ( i , j , n-1 >  +  ABt  ,  n=l  (Ill.Sa) 

T  1  I  T 

T1  (  i  ,  j  ,  n>  =  -  C-T1  (  i,  j  ,n-2)  +  AT^i.j.n-l)  +  2ABj]  ,  n>l  (Ill.Sb) 


where 


ABt  -  - -  1 2C T  ( i , j-1 , n-1 >  -  T  < i , j , n-1 ) ]/Ay2 


♦  CT1 ( i+1 , j , n-1 )  -  2T1 ( i , j , n-1 )  +  TJ( i-1 , j , n-1 ) /Ax2) . ( III . S) 


The  dimensionless  forms  were  given  in  Equations  (3.40a,b). 
On  face  AC  (Figure  I I I. 2) 


Ay.+Ay2  TII<i-l,j)  -  T1 1 < i , J > 
Ow-»p  =  kn(  2  ‘  Av 


0S->P  =  k^CAx/Z) 


On-»p  =  ^jj(Ax/2) 


TIr<i,j-l)  -  T11*!,;)) 


TII<i,  j+1)  -  T^Ci.J) 


(III. 7) 


(III. 8) 


(III. 9) 


UP  =  P  c  [Ax-(Ay1+Ay2>/41 


T1 1  (  i , j , n ) -T1 1 ( i , J  ,  n-1  > 


,  n=l  (III. 10a) 


up  =  p  c  IAx- (Ay^Ay^MI 


3T1 1 ( i , j , n )-4TT 1 ( i , j , n-1 >+Tr 1 < i , j , n-2  > 


n>l  (III. 10b) 


Up  =  09Om  gives 


T11 ( 1 , j ,n)  =  TII(i, j,n-l )  +  ABZ  , 


n=l  (III. 11a) 


T1 1 ( i , j , n)  =  -  [-T1I(i,j,n-2)+4TII(i,j,n-l >+2AB2 J ,  n>l  (III. lib) 


where 


w 

A 


ItM 


Hi  r ,,.11.  _ 


II  II 


ftT  (i-l.j.n-1)  -  T  (i,.1,n-l)]/Ax2 


+  ^  -  T  ( i, j , n-1 ) ]/(Ayt2  +  Ay1Ay2) 


II,,  , 


+  tT  (i*j+l.n-l)  -  T  (i, j,n-l)]/(Ay1Ayz  +  Ay22)} 


(III. 12) 


The  dimensionless  forms  were  given  in  Equations  <3.42a,b>. 
fti  face  BD  (Figure  I I I. 3) 


<w  •  k  I— 

KP  li  2  Ax 


Qs+p  ~  k  (Ax/2) 


TIr<  i  .  j~l )  -  T^d.J) 


(III. 13) 


(III. 14) 


Qn-*p  =  k^(Ax/2) 


T  ( i , j+1 )  -  T 1 1 ( i , j ) 


(III. IS) 


Up  =  PnciltAx.(Ay1^y2,/4]  J.h-1)  ^  n=1  (m>16a) 


UP  a  PI1cIjIAx'<Ay1 +Ay 2 ) /4 3 


UP  =  °av™  gives 


21 _ (1, J,n)-4T  (i, j,n-l)+Trr(i, j,n-2) 

2At  ~ 

n>l  (III. 16b) 


T  (i.j.n)  =  Tn(l,j,n-l)  +  AB3  , 


n=l  (111.17a) 


T  (i.J.n)  =  3  C-T  (i, J,n-2)+4TII( i, j,n-l)+2AB3I  ,  n>l  (III. 17b) 


T1 1 ( i , j , n)  =  ^  t-T1 1 ( i , j , n-2)+4T1 1 ( i , j , n-1 J-^AB^l  .  n>l  (III. 23b) 
where 

ABU  =  - —  l[Tn(i-l,  J,n-1>  -  2T1 1  (  i  ,  j ,  n-1 )  +  T1 1  (  i+1 ,  j ,  n-1 )  3/Ax2 

PIICII 

+  2CTri<i, j+l,n-l)  -  T11 ( i , j , n-1 >3}  .  (III. 24) 


The  dimensionless  forms  were  given  in  Equations  (3.46a,b). 

III. 2  Energy  Balance  at  the  Corner  of  the  Cavity  (Figure  3.3) 

If  one  has  a  two-dimensional  configuration  that  looks  like  the 

cross  section  of  two  walls  meeting  at  a  corner,  the  nodal  point  at  the 

corner  so  formed  is  called  a  reentrant  corner. 

Corner  A.  The  node  P  in  Figure  I II. 5  is  located  at  a  reentrant 

corner.  The  nodes  N,  R,  S,  and  W  are  called  exterior  corner  nodes. 

Observing  that  the  areas  available  for  the  flow  of  heat  from  N  to  P 

Ax  Ay 

and  from  R  to  P  are  *1  and  —  -1,  respectively.  The  heat  fluxes 
toward  the  corner  point  P(i,j>  are: 


kj+kj  x  T< 1-1 , J  >  -  T( i , J  > 

"  <_ 2~i->  <Ay> - Ax - 


<  III .  255 ) 


T( i , j+1 )  -  T(i,j) 

Qn->p  =  kjj<Ax/2)  ^ 


(III. 26) 


pi 


T< i+1 , j  )  -  T( i , j  > 

Or-»p  -  ki(Ay/2) 


(III. 27) 


Qc-)p  —  k  ( Ax ) 


T(i, j-1 >  -  T< i , J ) 


(III. 28) 


81 


**9 


-iso- 


*88 


l  e.  A  « 


>  .S' 


and  the  rate  of  change  of  internal  energy  is 


J 


ft 


I 


1:5 


-  ,  T( i , 1 , n)-T< i ,  i ,n-l > 

Dp  =  +  piicii/4)  <Ax‘Ay>  - ^ -  ,  n=l  (III. 29a) 


•  ,  3T( i , j , n)-4T( i , j ,n-l >+T< i , j ,n-2> 

DP  «  <PlCi/2  ♦  PllCn/4)  (Ax-Ay) - — - * - 


n>l  (III. 29b) 


From  conservation  of  energy,  Up  =  0sum,  one  obtains 


T(i,j,n)  =  T(i,j,n-1)  +  ABS  , 


n=l  (III. 30a) 


T(i, j,n)  =  -  I -T ( i , j , n-2 )  ♦  4T(i,j,n-l)  +  2AB5]  ,  n>l  (III. 30b) 


where 


Jjv  /  t-. 


I  I  'll  II 


+  k^tm,  j+l,n-l)-T(i,  j,n-l)l/Ay2  +  k  [T(i+1,  j,n-l) 


-  T(i, J,n-l)]/Ax2+2ki[T(i, j-l ,n-l >-T( i , J , n-1 ) I/Ay2} .  (III. 31) 


Corner  B  (Figure  I I I. 6) 


Qy->P  =  K  ^  ( Ay  /  2 ) 


T(i-l.j)  -  T( i , j ) 


(III. 32) 


Qn-»p  =  kjj(Ax/2) 


T(i,J+l>  -  T( i , j  > 


(III. 33) 


ki*kj j 

Qr-*p  =  <2*  <Ay) 


T<i-H,j)  -  T(i,J) 
Ax 


(III. 34) 


DTTEBFACE 


Figure  III. 6  Energy  balance  at  comer  B 


-J*  -Ji-jr-jrnjr 


k  .A*.  n^-1’  -T<1^. 

1  Ay 


(III.3S) 


0-  “  ‘V/2  *  A^/AKAx-Ay.  ^  ( 


III. 36a) 


Up  "  <PiCi/2  +  PI1cII/^)(Ax-Ay) 


3T<i, j,n)-4T(i, j,n-l>+T(j. j,n-2) 
2At 

n>l  (III. 36b) 


UP  =  Q«um  gives 


T<i> j»n)  =  T<i, j,n-l)  +  AB,  , 


n=l  (III. 37a) 


Td.j.n)  =  3  t-T(l,J,n-2>  +  4T(i,J,n-l)  ♦  2ABfel  ,  „>i  (111.37b) 


where 


AB  —  uv/ 

6  "  <pici/2+piicii/4)  {kiCT<1*lt ~  T(i, j,n-l)l/Axz 


+  kI1CT<itJ+l.n-l)  -  T(i, J,n-1)]/Ay2  +  (kj+krI >IT< i+1 , j , n-1 ) 
-  T< i » J t n-1 ) ]/Axz+2kj tT( i , J-l , n-1 )-T< i , J , n-i , ]/Ay2}  .(III. 38) 


Corner  C  (Figure  1 1 1. 7) 


Qw-*p  =  kn<Ay 


>  T  (i-i, j)  ~  TI1(j, j) 
Ax 


<111.39; 


SWISS' 


(III. 40) 


,11 


Qn-»p  =  kjj(Ax) 


T“<i,J+l>  ~  Tri(i,j> 
Ay 


.rr. 


Qr-»p  =  ^ll<-by/2) 


0g-»p  -  k^fAx/Z) 


T^(i4j,j)  -  Tn<i,j) 
Ax 


.XI 


T**<  i  i  j~l  >  -  Tn(i,j) 
Ay 


(III. 41) 


(III. 42) 


and 


Up  =  PHCli(4  Ax>Ay> 


3  TlI(i» J.n)  -  TI1(i,i,n-l) 


At 


n=l  (III. 43a) 


Up  =  piicii(4  Ax  Ay) 
Up  =  Osum  gives 


3Tir<l, jtn)-4T11 (1, j,n-l )»Tri(i, j,n-2) 
2At  ~ 


n>l  (III. 43b) 


■  ii. , 


T* * ( i  >  J » n)  =  Tn(i,  j.n-l)  +  AB-,  , 


n=l  (III. 44a) 


T  d.j.n)  =  3  t-T1I(i,j,n-2)+4TI1<i,j,n-l)+2AB.,]  ,  n>l  (III. 44b) 


where 


2AtkI}  TT 

AB?  =  ^PllCii  fC2T  * 1-1 >  J  »n_1 >“3Tr 1 ( i , J ,  n-1 l+T11 ( i+i f j, n-1 ) I/Ax2 


.II 


+  IT“(i, j-ifn-l)  -  3TXI<l,j,n-l)  +  2T 


.II 


( 1 , j+l,n-l ) ]/Ay2}  . 


(III. 45) 


Corner  D  (Figure  I I I. 8) 


Qw>p  -  kI1<Ay/2> 


T1 1  <  i-1 ,  j  >  ~  Tn(i,j> 
Ax 


II 


T1 1  ( i , j+1 )  -  T1 1 ( i , j  > 
Ay 


Qr-*p  -  kjj<Ay> 


T1I(i*l,J)  -  T11 1 i , j  > 
Ax 


Qs+p  ~  kij(Ax/2> 


Tir(i,j-1)  -  T1 1 ( i , j ) 
Ay 


and 


°P  =  PIICII(4  Ax'Ay> 


T 1 1 ( i , j , n)  -  T11 ( i , j , n-1 ) 
At 


n-1 


Up  -  p  c  (-  Ax. Ay) 
p  Kil  II  4 


3T1 1 ( i , j , n ) -4T 1 1 ( i , j , n- 1 >  +T 1 1 ( i , j , n- 


2At 


n>l 


UP  =  08um  8ives 


TII<1, j,n)  *  TII(i,j,n-l)  +  AB8  ,  n=l 

Tn(i,j,n)  =  3  t-T11^,  j,n-2)+4TI]I<i,j,n-l)+2AB8]  ,  n>l 


where 


kr  t 

AB« =  3r~r" 

ii 


f  CT11 < i-1 , J ,n-l )-3TJ  J< i , j  ,n-l )+2TTI( i+1, j , 


< III. 46) 

(111. 47) 

(111. 48) 

(111. 49) 

(111.50a) 

•2) 

(111. 50b) 

(111.51a) 

(III. Sib) 

n-1 ) I/Ax2 


+  CT1 1 ( i , 1-1 , n-1 )  -  3T1 1 ( 1 , j , n-1 )  +  2T11 ( i , j+1 , n-1 ) I /Ay2] . ( II I . 52 ) 


The  dimensionless  forms  for  the  corner  points  A,  B,  C,  and  D  were  given 
in  Equations  <3.48a,b,  3.S0a,b,  3.S2a,b,  3.S4a,b>,  respectively. 


* 

■  i" 


i* 


III. 3  Energy  Balance  on  the  Surface  Boundary  (Figure  III. 9) 

For  the  explicit  scheme,  the  time  step  is  limited  by  the  stability 
criterion.  As  a  result,  the  moving  asperity  at  some  time  may  not  be 
right  above  the  grid  points.  To  alleviate  this  situation,  one  can 
also  use  the  energy  balance  method  to  describe  the  surface  boundary 
condition. 

The  heat  fluxes  toward  P(i,j)  from  material  points  N,  R,  and  W  are 
Tr(i-l,l)  -  T1 ( i  ,1 > 

Qw->p  =  k^Ay/2)  - — -  ,  (III. S3) 

TI(i,2)  -  TI(i,l) 

Qn-»p  =  kj(Ax  ’  Ay  (III.S4) 

T1  ( i+1 , 1 )  -  T1 ( i , 1 > 

Qr->p  =  k^Ay/2)  - - - .  (III.  SS > 

The  exterior  heat  which  is  conducted  into  the  neighborhood  surface  of 
the  boundary  point  P(i,J),  which  is  under  the  asperity,  is 

Qext.  =  q<0.5Ax  +  h)/unit  thickness  ,  (III. 56) 


4' 


f 


where  h  is  less  than  Ax/2.  The  formulation  thus  takes  care  of  all 
cases  when  the  asperity  end  points  do  not  fall  on  the  grid  point.  The 
rate  of  change  of  the  internal  energy  Up  in  the  interval  At  at  P(i,j) 
is 


Up  = 


PjC^Ax-Ay/2) 


TI(i,J,n)  ~  T 1 < i , j , n-1 ) 
At 


n=l  (III. 57a) 
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UP  =  Pjcj <Ax- Ay/2) 


3T1 < i , j , n)  -  4TI(i,j,n-l)  +  TI(i,J,n-2) 


2At 


n>l  (111.57b) 


From  conservation  of  energy,  Up  =  08um,  one  obtains 


T  r  2q(Ax/2+h)At 

T  (i,l,n>  =  T  <i,l,n-l)  +  —  — — -  +  ABU  , 

p  c  AxAy  y 


n=l  (111.58a) 


I  I 


T  < i , 1 , n) 


2q< Ax/2+h)At 


5  C_T  +  4T  ‘i.j.n-D  *  p  c  ^ 


+  2ABy], 


I  I 


n>l  (111.58b) 


where 


Atk, 


ABq  =  - 1  ([TJ(i-l,l,n-l )  -  2TI(i,l,n-l)  +  TI(i+l,l,n-l)l/Ax2 

’  p  c 

i 


+  2[TI(i,2,n-l)  -  TI(i,l,n-l)]/Ay2j  . 


(1II.S9) 


The  dimensionless  forms  for  the  surface  boundary  were  given  in 
Equations  <3.S6a,b). 


-161- 


APPENDIX  IV 


THE  PROGRAMS  TO  COMPUTE  THE  TEMPERATURE  AMD  THE  STRESS 


FIELDS  SOLITTIONS 


C  MAIN  PROGRAM 

C  TEMPERATURE  FIELD  OF  A  LAYERED  MEDIUM  WITH  A  CAVITY,  THE  TOP 
C  EDGE  OF  THE  CAVITY  IS  AT  THE  INTERFACE 
C 

IMPLICIT  REAL* 8  <A-H,0-Z> 

C 

C  T,TT,TTT  =  TEMPERATURE  IN  THE  CURRENT,  PREVIOUS  ONE,  AND 
C  PREVIOUS  TWO  TIME  STEPS 

C  FX  &  FY  =  HEAT  FLUX  IN  U  Y  DIRECTION,  RESPECTIVELY 
C  0  =  SURFACE  HEAT  INPUT 

C  X  &  Y  =  COORDINATES  IN  PHYSICAL  PLANE 

C  XS,YE,XSS, YEE _  =  THE  DERIVATIVES  OF  THE  COORDINATES  IN 

C  PHYSICAL  PLANE  WITH  RESPECT  TO  THE  COORDINATES  IN  COMPUTATIONAL. 

C  PLANE 

C 

DIMENSION  T( 14S, 33) ,TT( 145, 33  > , TTT(145, 33) , FX(145, 33) ,FY(145, 33) , 
AQ< 14S ) , X( 14S , 33 ) , Y< 145 , 33 ) , XS( 14S , 33 ) , YE( 14S , 33 ) 

C 

C  DA,DG,SIG,TAU  =  THE  COEFFICIENTS  DEFINED  IN  TEMPERATURE  EQUATION 

C 

DIMENSION  XSS<14S,33) , YEE <145, 33) ,DA(14S,33) ,DG(1 45,33 ) , 

ASIG< 14S , 33  > ,TAU( 14S , 33 ) 

C 

C  TS,FXS,FYS  =  TEMPERATURE  AND  ITS  GRADIENTS  IN  THE  CORRESPONDING 
C  POINTS  IN  STRESS  FIELD 
C 

DIMENSION  TS<  67 , 35 ) , FXS<  67 , 35 ) , FYS<  67,35) 

COMMON  /AS4/  RL1 , RL2 , RLI , RLC , RMU1 , RMU2 , RMUI , RMUC , EX1 , EX2 , 

AEX I , EXC , HAK , RHC 

COMMON  /ASS/  DV,DL,DIFF1,DIFF2, CONDI, COND2, BETA, R1,R2, 

AALPHS , DX , DY 1 , DY2, DT, Rll ,R12, R21 ,R22, A1 , A2, A3, A4, AS, A, AA, A6, A7 , A8 
COMMON  /A56/  DTR1 , DTR2 , Ml , M2 , N1 , N2 , MM1 , MM2 , NN1 , NN2 , NA1 , NA2 , 

AMA1 , MA2 , K1 , K2 , K3 , NTE , M122 , M222 , M112 , M212 , N121 
COMMON  /AS7/  ID2, ID21 , ID3, ID31 , 12, J2, 13, J3, 14, J4 
C 

C  II  &  J1  =  TOTAL  GRID  POINTS  OF  THE  TEMPERATURE  FIELD  IN  X  &  Y 
C  DIRECTION,  RESPECTIVELY 

C 

11=145 
Jl  =  33 
C 

C  LI1  &  KJ1  =  TOTAL  GRID  POINTS  OF  THE  STRESS  FIELD  IN  X  8.  Y 
C  DIRECTION,  RESPECTIVELY 

C 

LI1=67 
KJ1=35 
12=11-1 
13=11-2 
J2=J1-1 
J3=Jl-2 
14=11-3 
J4= Jl-3 
C 

C  Ml  *■  M2  =  X  COORDINATES  OF  THE  CAVITY  CORNERS 
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Ml  =  33 
M2=t>4 


N1  &  N2  =  Y  COORDINATES  OF  THE  CAVITY  CORNERS 


Nl=14 

N2=24 

MM1=M1-1 

MM2=M2+1 

NN1=N1-1 

NN2=N2+1 

NA1=N1+1 

NA2=N2-1 

MA1=M1+1 

MA2=M2-1 

M122=Ml+2 

M222=M2+2 

M112=Ml-2 

M212=M2-2 

N121=N1+1 


K1  =  LAYERED  THICKNESS 


K1=N1 

K2=K1-1 

K3-K1+1 

ID2=8 

ID21=ID2+1 

ID3=140 

ID31=ID3-1 


NTE  =  FINAL  TIME  STEP 


NTE=121 


DV  =  TRAVERSE  SPEED  OF  ASPERITY 


DV=6.D2 


DL  =  ASPERITY  WIDTH 


DL=1 . D-2 


COND  =  THERMAL  CONDUCTIVITY 


CONDI =1. 21 3D0 

COND2=l . 213D0 

COND I = ( CONDI ♦COND2 ) / 2 . DO 

CONDC  =  <  2 . DO*COND1+COND2  >/3 . DO 


DIFF  =  THERMAL  DIFFUSIVITY 


DIFF1=4. 29D-3 


DIFF2=4. 29D-3 
DIFFI=(DIFF1+DIFF2)/2.D0 
DIFFC=(2.DO*DIFFl«-DIFF2)/3 . DO 
C 

C  RNU  =  POISSON'S  RATIO 
C 

RNU1=0. 285D0 
RNU2=0. 28SD0 
RNUI=r(RNUl+RNU2)/2.D0 
RNUC=(2.D0*RNU1+RNU2)/3.D0 
C 

C  E  =  YOUNG’S  MODULUS 

C 

E1=3.6D7 
E2=3.6D7 
El- <E1+E2)/2.D0 
EC= ( 2 . DO*El+E2 ) /3 . DO 
C 

C  RHO  =  MASS  DENSITY 

C 

RH01=9. 31D-3 
RH02=9. 31D-3 
RHOI=< RHOl +RH02 ) / 2 . DO 
RHOC=( 2 . DO*RH01+RH02 ) /3 . DO 
C 

C  EX  =  THE  COEFFICIENT  OF  THERMAL  EXPANSION 
C 

EX1=6. 29D-6 
EX2=6. 29D-6 
EX I = ( EX1+EX2 ) /2 . DO 
EXC=<  2 . D0*EX1+EX2 ) /3 . DO 
C 

C  RMUF  =  COULOMB  COEFFICIENT  OF  FRICTION 
C 

RMUF=0. 5D0 
C 

C  RMU  &  RL  =  LAME  CONSTANTS 
C 

RMU1=E1/(2.D0*<1.D0+RNU1 ) ) 

RMU2=E2/ « 2 . DO*  < 1 . D0+RNU2 ) ) 
RMUI=EI/(2.D0*<1 .DO+RNUI >  > 

RMUC=EC/ ( 2 . DO* ( 1 . DO+RNUC ) ) 

RL1 =2 . D0*RMU1*RNU1 / ( 1 . DO-2 . D0*RNU1 > 

RL2=2 . D0*RMU2*RNU2/ < 1 . DO-2 . D0*RNU2  > 
RLI=2.D0*RMUI*RNUI/< 1 .D0-2.D0*RNUI ) 

RLC=2. DO*RMUC*RNUC/ < 1 . DO-2 . D0*RNUC ) 
C2=DSQRT ( RMU2* 1 . 2D1 /RH02  > 

HAK =RMUF • DV * DL/COND 1 

RHC=RH02*C2**2 

BETA=COND2/CONDl 

R1=DV*DL/DIFF1 

R2=DV*DL/DIFF2 

ALPHS=DIFF2/DIFF1 

DX=0 . 2D-1 


fib 

W- 

m 

m.* 


DY2=0. 2D-1 
C 

C  DT  =  TIME  STEP 
C 

DT=0.1D-1 

Ril=DT/<Rl*DX*DX) 

R12=DT/(R1*DY1*DY1 ) 

R21=DT/<R2*DX*DX> 

R22=DT/(R2*DY2*DY2) 

A1 =1 . DO+BETA/ALPHS 

A2=1.D0+BETA 

A3=A2/A1 

A4=l . DO-2 . DO*A3*Rll-2 . D0*A3*R12 
AS=0. SDO+O. 2SD0*BETA/ALPHS 
A=1 . DO-2 . D0*R21-2 . D0*R22 
AA=1 . DO-2 . D0*Rll-2 . D0*R12 

A6=0 . SD0*A2*Rll+0 . 5D0* BETA* R1 2+0 . SD0*R11+R12 
A7=R1+BETA*R2 
A8=A2/A7 
DTR1=DT/R1 
DTR2=DT/R2 
C 

CALL  XYLCTCX, Y,XS, YE,XSS, YEE,DA,DG,SIG,TAU, II , J1 > 

C 

HC=1 

NT=1 

100  WT=NT+1 

TIME=(NT-1)*DT 

C 

CALL  QIN(DX,DT,NT, II ,Q, II , JJ) 

C 

CALL  TEMP ( T , TT , TTT , 0 , FX , FY , X , Y , XS , YE , XSS , YEE , DA , DG , SIG , TAU , 
All , JJ.MT, II , J1 > 

C 

NC=NC+1 

IF(NT  .HE.  NTE)  GO  TO  S 
C 

CALL  MAP(TS,FXS,FYS,T,FX,FY, II , J1 ,LI1 ,KJ1 ) 

C 

DO  31  J=NN1,NA1 
WRITE(6,*)  J 

WRITE< 6, * )  (T< I , J) , I=MM1 , MAI ) 

WRITE<6, * )  ( FX ( I , J > , I =MM1 , MAI > 

31  WRITE(6, * )  ( FY ( I , J ) , I =MM1 , MAI ) 

DO  32  J=NNi ,NA1 

WRITE< 6, * )  J 

WRITE(6, *)  (T( I, J> , I=MA2,MM2) 

WRITE<6,*)  <FX< I , J) , I=MA2,MM2) 

32  WRITE(6, * )  ( FY ( I , J ) , I =MA2 , MM2 > 

DO  33  J=NA2,NN2 

WRITE(6, • )  J 

WRITE< 6, * )  ( T ( I , J ) , I =MM1 , MAI ) 

WRITE<6, * )  ( FX ( I , J ) , I =MM1 , MAI ) 


t 
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33  WRITE<6, * )  <  FY  < I , J  > , I =MMi , MAI ) 
DO  34  J=NA2,NN2 

WRITE< 6, • )  J 

WRITE( 6, * )  ( T ( I , J > , I =MA2 , MM2 ) 
WRITE<  6, • )  < FX < I , J ) , I =MA2 , MM2 ) 

34  WRITE<6,*>  ( FY ( I , J ) , I =MA2 , MM2 ) 
WRITE(6,101) 

101  FORMAT  < / , SX , ' TS ’ , / ) 

DO  20  J=1,KJ1 

20  WRITE<6,90>  (TS(I,J) ,1=1, LIl) 
WRITE(6,102) 

102  FORMAT  < / , 5X , * FXS '  ,  / ) 

DO  21  J=1,KJ1 

21  WRITE<6, 90)  (FXS< I , J ) , 1=1 , LIl ) 
WRITE<6,103) 

103  FORMAT ( / , SX , ' FYS ' , / ) 

DO  22  J=1 , KJ1 

22  WRITE<6,90>  (FYS( I , J ) , 1=1, LIl ) 
90  FORMAT  <  S ( IX , D14 . 7 ) ) 

S  IF(NT  .LT.  NTE)  GO  TO  100 
STOP 


C  THIS  SUBROUTINE  INPUTS  COORDINATES  X  &  Y  AND  CALCULATES  THE 

C  DERIVATIVES  OF  THE  COORDINATES 

C 

SUBROUTINE  XYLCT<X,Y,XS,YE,XSS,YEE,DA,DG,SIG,TAU, II, Jl) 
IMPLICIT  REAL *8  (A-H.O-Z) 

C 

C  X  &  Y  =  COORDINATES  IN  PHYSICAL  PLANE 
C  XS , YE , XSS , YEE  =  THE  DERIVATIVES  OF  THE  COORDINATES  IN 

C  PHYSICAL  PLANE  WITH  RESPECT  TO  THE  COORDINATES  IN 

C  COMPUTATIONAL  PLANE 

C  DA.DG.SIG.TAU  =  COEFFICIENTS  DEFINED  IN  TEMPERATURE  EQUATION 

C 

DIMENSION  X<I1,J1),Y(I1,J1),XS(I1,J1),YE(I1,J1) ,XSS( II , Jl ) 
DIMENSION  YE£< II , Jl ) , DA( II , Jl ) ,DG( II , Jl ) ,SIG< II , Jl ) ,TAU< II , Jl > 
12=11-1 
J2=J1-1 
13=11-2 
J3=Jl-2 
C 

DO  90  J=1,J1 
DO  31  1=1,11 
X(I,J)=0.D0 
Y< I,J)=0.D0 
XS(I,J)=0.D0 
YE< I,J)=0.D0 
XSS< I , J)=0.D0 
YEE( I , J)=0.D0 
DA( I , J )=0.D0 
DG< I, J)=0.D0 
SIG( I , J)=0.D0 
TAU( I,J)=0.D0 
91  CONTINUE 
90  CONTINUE 
C 

DO  1  J=1,J1 
C 

DO  5  1=1,4 

5  X( I , J)=( 1-1 )*0 . SDO 
C 

DO  6  1=5,9 

6  X( I , J ) =1 . 5D0+< 1-4) *0. IDO 
C 

DO  7  1=10,134 

7  X( I,J)=2.D0+( I-9>*O.O2D0 
C 

DO  8  1=135,138 

8  X( I , J  >  =4. SD0+( 1-134 )*0 . OSDO 
C 

DO  9  1=139,145 

9  X( I, J)=4.7D0+< 1-138 )*0. IDO 
C 

1  CONTINUE 
C 

DO  3  1=1,11 


t  t  wi 
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Y  ( 1 , 1 >  =0 .  DO 

Y(  1 ,2)=0.00SD0 
Y< I,3>=0.01D0 
Y ( 1 , 4)=0. 016D0 
Y(  I , S)=0. 023D0 
Y< I,6)=0.031D0 
Y( I , 7 )=0. 04D0 
Y< 1 , 8 )=0. 0SD0 

Y  < 1 , 9) =0. 059D0 
Y( I , 10)=0. 067D0 
Y( I,11)=0.07SD0 
Y< I , 12)=0. 082D0 
Y< I,13>=0.088D0 
Y( I ,14>=0. 094D0 
Y( I,1S)=0.1D0 
Y< I , 16)=0. 107D0 
Y< I , 17 ) =0. 11SD0 
Y< I, 18>=0.125D0 
Y( I,19)=0.137D0 
Y( I,20)=0.1S1D0 
Y( I , 21 )=0. 167D0 
Y( I , 22)=0. 185D0 
Y< I,23)=0.205D0 
Y( I , 24)=0. 22SD0 
Y< I , 2S) =0. 24SD0 
Y< I , 26>=0. 27D0 
Y(  I , 27 ) =0 . 3D0 
Y( 1 , 28)=0. 34D0 
Y< I,29)=0.4D0 
Y(  I , 30)=0. 48D0 
Y( I , 31 )=0. S8D0 
Y( I , 32 ) =0. 78D0 
Y( I , 33 ) =1 . 08D0 

3  CONTINUE 
C 

DO  17  J-2, J2 
DO  18  1=2,12 

XS(I,J)=(X(I+1,J)-X<I-1,J) )/2.D0 
YE(I,J>=<Y(I,J+1>-Y(I,J-1))/2.D0 
XSS( I , J )=X< 1-1 , J  >-2. D0*X< I , J )+X( 1+1 , J ) 

YEE( I , J) =Y( I , J-l >-2 . D0*Y ( I , J )+Y( I ,  J+l ) 

18  CONTINUE 
17  CONTINUE 

C 

1=1 

DO  19  J=2,J2 

XS< I , J)=(-3.D0*X( I , J)+4.D0*X< 1+1 ,J)-X( 1+2, J ) )/2.D0 
YE< I , J ) =< Y< I , J+l )-Y ( I , J-l > )/2. DO 

19  CONTINUE 
C 

J=1 

DO  20  1=2,12 

XS< I , J )  =  <  X< 1+1 , J  >-X( 1-1 , J ) )/2. DO 

YE(I,J)  =  (-3.D0*Y( I,J>+4.D0*Y< I, J+l )-Y( I , J+2 >  >/2 . DO 


CONTINUE 


1  =  11 

DO  21  J=2,J2 

XS(I,J)  =  (X(I-2,J>-4.D0*X(I-1,J>+3.D0*X(I,J>  >/2.D0 

YE( I , J)=< Y< I ,  J+l )-Y< I , J-l ) )/2. DO 

CONTINUE 


J=J1 

DO  22  1=2,12 

XS<I,J)  =  <X<I+1,J)-X(I-1,J>  >/2.D0 

YE< I , J  >  =  ( Y( I , J-2 )-4. D0*Y< I , J-i )+3.D0*Y ( I , J) )/2 . DO 

CONTINUE 

XS( 1 , 1 )=XS( 1 ,2) 

XS<1,J1)=XS<1,2) 

XS( 11,1 )=XS< 11,2) 

XS( II, J1)=XS( 11,2) 

YE< 1 , 1 )=YE<  2, 1 > 

YE(1,J1)=YE(2,J1) 

YE< 11,1 )=YE(2,1 > 

YE< II , J1 )=YE( 2, J1 > 

DO  2S  J=1,J1 

DO  28  1=1,11 

DA< I , J)=l .DO/XS< I , J  >**2 

DG( I , J) =1 .D0/YE< I , J>**2 

SIG< I , J)=-YEE( I , J)/YE< I , J)**3 

TAU( I,J)=-XSS( I,J)/XS< I,J)**3 

IF<DABS<SIG<  I, J))  .LT.  l.D-10)  SIG( I , J)=0.D0 

IF<DABS<TAU<  I, J))  .LT.  l.D-10)  TAU< I , J)=0.D0 

CONTINUE 

CONTINUE 

RETURN 

END 


c 

c 


THIS  SUBROUTINE  INPUTS  THE  SURFACE  B.C. 


as 


SUBROUTINE  QIN(DX,DT,NT, Ii,Q, II , JJ) 
IMPLICIT  REAL *8  (A-H.O-Z) 

0  =  SURFACE  HEAT  INPUT 

DIMENSION  Q< II ) 

NR2=DX/DT,NRl=NR2/2 

NR1=1 

NR2=2 


DO  24  1=1,11 
Q(  I )=0.D0 
CONTINUE 


m 


8 

& 


IF(NT  .GT.  2)  GO  TO  21 

NC1=0 

NC2=0 

NC3=0 

NTT=1 

ND1=NR1 

ND2=NR2 

ND3=0 

NNT=1 

NC1=NC1+1 

NC2=NC2+1 

IF (NCI  .EQ.  NR1 )  GO  TO  1 

GO  TO  4 

NC3=NC3+1 

IF<NC3  .EQ.  1>  GO  TO  6 

NC3=0 

GO  TO  7 

NTT=NTT+1 

NC1=0 

II=74-NTT+1 

IF(NC2  .GE.  NR1)  GO  TO  2 

PB=O.SDO 

GO  TO  5 

IF(NC2  .LT.  NR2)  GO  TO  3 

PB=0. SDO 

NC2=0 

GO  TO  S 

PB=0.D0 

0<  II) =PB+DT*NC1/DX 

ND1=ND1-1 

ND2=ND2-1 

IF(ND1  .EQ.  0)  GO  TO  11 

GO  TO  12 

ND3=ND3+1 

IF(ND3  .EQ.  1)  GO  TO  13 
ND3=0 


>3 


GO  TO  14 
NNT=NNT+1 
ND1=NR1 
JJ=124-NNT+1 

IFCND2  .LE.  NR1 )  GO  TO  IS 

PF=0.D0 

GO  TO  16 

IF<ND2  .GT.  0)  GO  TO  17 

ND2=NR2 

PF=0.D0 

GO  TO  16 

PF=0. SDO 

Q<JJ)=PF+DT*ND1/DX 

111=11+1 

JJ1=JJ-1 

DO  23  I=II1,JJ1 

0<I)=1.DO 

CONTINUE 

RETURN 

END 


c 

c 

c 


THIS  SUBROUTINE  SOLVED  A  LAYERED  MEDIUM  WITH  A  CAVITY,  THE  TOP 
EDGE  OF  THE  CAVITY  IS  AT  THE  INTERFACE 


SUBROUTINE  TEMP < T , TT , ITT , Q , FX , FY , X , Y , XS , YE , XSS , YEE , DA , DG , 
ASIG.TAU, II,JJ,NT,I1,J1> 

IMPLICIT  REAL *8  (A-H.O-Z) 

DIMENSION  T(Il.Jl) ,TT<I1,J1),TTT(I1,J1) ,FX(I1,J1>, 

AFY  <I1,J1),Q(I1),X<I1,J1>,Y(I1,J1),XS(I1,J1),YE<I1,J1) 

DIMENSION  XSS< II , J1 > , YEE( I1,J1),DA(I1,J1),DG( II, Jl) , 

ASIG( II , J1 ) ,TAU< II , J1 ) 

COMMON  /AS4/  RLi , RL2 , RLI , RLC, RMU1 , RMU2 , RMUI , RMUC , EXl , EX2 , 

AEX I , EXC , HAK , RHC 

COMMON  /ASS/  DV,DL,DIFF1,DIFF2, CONDI, COND2, BETA, R1,R2, 

AALPHS, DX, DY1 , DY2, DT, Rll ,R12, R21 , R22 , A1 , A2, A3, A4, AS, A, AA, A6, A7 , A8 
COMMON  /AS6/  DTR1 , DTR2 , Ml , M2 , N1 , N2 , MM1 , MM2 , NN1 , NN2 , NA1 , NA2 , 

AMA1 , MA2 , K1 , K2 , K3 , NTE , M122 , M222 , M112 , M212 , N121 
COMMON  /AS7/  ID2, ID21 , ID3, ID31 , 12, J2, 13, J3, 14, J4 
IF(NT  .NE.  2)  GO  TO  100 
DO  1  J=1,J1 
DO  2  1=1,11 
T(I,J)=O.DO 
TT( I,J)=0.D0 
TTT< I , J )=0.D0 

2  CONTINUE 
1  CONTINUE 

C 

C  COMPUTE  TEMPERATURE  OF  THE  COATING  LAYER 

C  DA,DG,SIG,TAU  =  THE  COEFFICIENTS  DEFINED  IN  THE  TEMPERATURE 

C  EQUATION 

C 

100  DO  3  1=2,12 
DO  4  J=2,K2 

AA1=DA< I,J)*(TT<I-1,J)-2.D0*TT(I , J )+TT< 1+1 , J ) ) 
AA2=DG(I,J)*(TT<I,J-1)-2.D0*TT( I , J >+TT< I , J+l > ) 

AA3=SIG<  I,J)*(TT<I,J+1  )_TT(  I ,  J-l )  >/2.D0 
AA4=TAU< I , J > * (TT( I+l , J )-ri( 1-1 , J ) >/2 . DO 
AAA=AA1+AA2+AA3-<AA4 
IF(NT  .GT.  2)  GO  TO  5 
C 

C  2  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 
C 

T( I , J ) =TT< I , J )+AAA*DT/Rl 
GO  TO  4 
C 

C  3  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 
C 

S  T< I , J)=(2.D0*AAA*DT/R1+4.D0*TT( I , J)-TTT( I , J) >/3.D0 
4  CONTINUE 

3  CONTINUE 
C 

C  COMPUTE  TEMPERATURE  OF  THE  SUBSTRATE 
C 

DO  b  1=2,12 
DO  7  J=K3,J2 
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AA1=DA( I , J)*(TT( 1-1 , J )-2 . DO*TT  < I , J  >+TT( I+l.J) > 
AA2=DG( I,J)*(TT(I,J-1 )-2. DO*TT < I , J )+TT< I , J+l >  > 
AA3=SIG( I,J)*<TT<I,J+1 )-TT  < I , J-l ) )/2. DO 
AA4=TAU( I,J)*<TT(I+1,J )-TT (I-1,J))/2.D0 
AAA=AA1+AA2+AA3+AA4 
IF(NT  .GT.  2)  GO  TO  8 


2  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 


T ( I , J)=TT( I , J )+AAA*DT/R2 


GO  TO  7 


3  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 


T( I , J ) = ( 2 . DO* AAA*DT/R2+4 . D0*TT ( I , J ) -TTT ( I , J ) ) /3 . DO 

CONTINUE 

CONTINUE 


DO  9  J=N1 ,N2 
DO  10  I=M1 , M2 
T<I,J)=0.D0 
CONTINUE 
CONTINUE 


COMPUTE  TEMPERATURE  OF  THE  CORNER  POINTS 


IF(NT  .GT.  2)  GO  TO  12 


2  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 


CORNER  A 


T(M1 ,N1 )=TT(M1 ,N1 )  +  <  0 . 5D0*A2*R11*TT(M1-1 ,N1 )- 
AA6*TT<  Ml , N1 ) +0 . 5D0*BETA*R12*TT( Ml , Nl+1 >+0 . SD0*R11*TT<  Ml+1 , N1 )  + 
AR12*TT  <  Ml , Nl-1 ) ) /AS 


CORNER  B 


T<M2,N1 >=TT<M2,N1 )+<0. 5D0*R11*TT<M2-1 ,N1 )-A6 
A*TT(M2,N1>+0.SD0*BETA*R12*TT<M2,N1+1>+0.SD0*A2*R11*TT(M2+1,N1) 
A+R12*TT  <  M2 , Nl-1 ) ) /AS 


CORNER  C 


T(M1 ,N2 )=A*TT<M1 , N2 )+4. DO* (R21*  <TT<M1-1 ,N2>  + 

AO . SD0*TT<  Ml+1 , N2 ) )+R22*  <  0 . SDO'TT ( Ml , N2-1 ) +TT( Ml , N2+1 ) ) ) /3 . DO 


CORNER  D 


T ( M2 , N2 ) =A*TT<  M2 . N2 ) +4 . DO*  <  R21 *  <  0 . 5D0*TT<  M2-1 , N2 ) 

A+TT<  M2+1 , N2 ) ) +R22*  <  0 . SD0*TT( M2 , N2-1 >+TT<  M2 , N2+1 ) ) ) /3 . DO 


GO  TO  14 
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3  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 


CORNER  A 

12  T<M1,N1)=(-TTT(M1,N1)+4.D0*TT(M1,N1 >+2.D0* 

A<0.SD0*A2‘R11*TT(M1-1,N1)-A6*TT(M1,N1)+0.SD0*BETA*R12*TT<M1, 
AN1+1 )+0. SD0*R11*TT<M1+1 ,N1 )+R12*TT(Ml ,N1-1 ) ) /AS)/ 3. DO 

CORNER  B 

T  <  M2 , N1 )  =  < -TTT ( M2 , N1 ) +4 . D0*TT  <  M2 , N1 )+2.D0* 
A(0.SD0*R11*TT<M2-1,N1)-A6*TT<M2,N1>+0.SD0*BETA*R12*TT(M2,N1+1 )+ 
AO . SDO* A2*R11*TT< M2+1 , N1 >+R12*TT ( M2 , Nl-1 ) ) /AS ) /3 . DO 

CORNER  C 

T(M1 , N2 )  =  <-TTT(Ml ,N2  >+4. D0*TT( Ml ,N2>+8. DO* 

A<  R21 *  <  TT <  Ml-1 , N2  > -1 . SDO *TT <  Ml , N2  >  +0 . SDO *TT <  Ml+1 , N2  >  >  +R22*  <  0 . SDO 
A*TT(M1,N2-1)-1.SD0*TT(M1,N2)+TT(M1,N2+1 ) ) )/3.D0)/3.D0 

CORNER  D 

T( M2 , N2 ) = ( -TTT ( M2 , N2 )+4 . D0*TT ( M2 , N2 ) +8 . DO* 

A(R21* ( 0. 5D0*TT(M2-1 ,N2  >-l . SD0*TT(M2 ,N2  >+TT(M2+l , N2 ) )+R22* (0. SDO* 
ATT(M2,N2-1)-1.SD0*TT(M2,N2)+TT(M2,N2+1> ) )/3.D0)/3.D0 

COMPUTE  TEMPERATURE  ON  THE  LEFT  &  RIGHT  HAND  EDGE  OF  THE  CAVITY 

14  DO  IS  J=NA1 ,NA2 
I=M1 

DDY1=Y( I , J)-Y< I , J-l ) 

DDY2=Y ( I , J+l )-Y( I , J ) 

YY1=DDY1**2+DDY1*DDY2 
YY  2 =DDY 1  *  DDY  2+DDY  2**2 

EE1=DT*< (TT( 1-1, J)-TT< I , J) >/DX**2+<TT< I , J-l )-TT( I , J ) )/YYl+ 

A(TT( I , J+l )-TT< I , J) )/YY2)/R2 
IF(DABS(EE1 )  .LT.  l.D-65)  EE1=0.D0 
I=M2 

EE2=DT*( (TT( 1+1 , J)-TT( I , J) )/DX**2+(TT( I , J-l )-TT( I , J) )/YYl+ 

A<TT< I , J+l )-TT< I , J ) )/YY2  >/R2 
IF(DABS<EE2)  .LT.  1.D-S5)  EE2-0.D0 
IF(NT  -GT.  2)  GO  TO  16 

2  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 


I  =M1 

T(I,J>=TT<I,J>+2.D0*EE1 
I -M2 

T( I , J ) =TT< I , J )+2 . D0*EF2 
GO  TO  15 

3  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 


16  I=M1 

T ( I , J)= (-TTT ( I , J)+4.D0*TT( I , J )+4 . DO’EEl )/3.D0 
I=M2 

T< I , J  >  =  <-TTT( I,J)+4.D0*TT(I,J )+4.D0*EE2)/3.D0 
15  CONTINUE 

COMPUTE  THE  TEMPERATURE  ON  THE  TOP  8.  BOTTOM  EDGE  OF  THE  CAVITY 

DO  17  I =MA1 , MA2 

IF (NT  .GT.  2)  GO  TO  18 

2  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 
J-Nl 

T  < I , J  >=AA*TT ( I , J )+Rll*  <TT( 1+1 , J )+TT( I-1,J))+2.D0*R12*TT(I,J-1 ) 
J=N2 

T( I , J ) =A*TT ( I , J )+R21* (TT< 1-1 , J  >+TT  < 1+1 , J) )  +  2.D0*R22*TT( I , J+l ) 
GO  TO  17 

3  POINTS  BACKWARD  DIFFERENCE  IN  TIME  DERIVATIVE 
18  J-Nl 

T ( I , J )  =  (-TTT ( I , J )+4.D0*TT( I , J )+4.D0* (Rll* ( 0. 5D0*TT  < 1+1 . J>- 
ATT ( I , J ) +0 . SD0*TT  < 1-1 , J) )+R12* (TT( I . J-l )-TT< I . J ) ) > >/3.D0 
J=N2 

T(I,  J)  =  (-TTT(I,  Jl+4.D0*TT(I,J>+4.DO*(R21*(O.5DO*nu-l .  J) 

A-TT ( I , J )+0 . SD0*TT< 1  +  1 , J  > )+R22* (TT ( I , J+l )-TT ( I . J  >  > >>/3.D0 

17  CONTINUE 

COMPUTE  THE  TEMPERATURE  AT  THE  INTERFACE 

DO  IS  IJ=1 , 2 

IF ( IJ  .EQ.  1)  GO  TO  20 

MN1=MM2 

MN2=I2 

GO  TO  21 

20  MNl=2 
MN2=MMl 

21  DO  22  I=MN1 . MN2 
J=K1 

DDX1=X( I,J>-X(I-1,J> 

DDX2--X<  1  +  1 ,  J)-X(  I  .J  ) 

XXI =DDX1 **2+DDXl *DDX? 

XX2rDDXl*DDX2+DDX2* *  2 
EE3=2.D0*A8‘DT* <TT<  1-1  ..J  »-TT  ■  i  X‘  • 

ATT ( I , J ) ) /XX2+2 . D0*DT*  <  TT  <  I  -1 
ABETA’DT*  (TT  < 1 . J  +  l  » - 71 1  I  .  ’  •'  V"  ** 

IF< NT  .GT.  2>  GO  TO  2.< 

2  POINTS  BACKWARD  *>  "• 


T(  I ,  J  >  TT (  1 ,  1 1 ♦(  l  :> 
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GO  TO  61 

57  EF121  =  <  3.D0*RLC+2. D0*RMUC>*1 . 2D1*EXC*HAK/RHC 

61  IF ( I  .LT.  MAI  .OR.  I  .GT.  MM2)  GO  TO  62 
IF < I  .GT.  MAI  .AND.  I  .LT.  MM2)  GO  TO  63 
IF( I  .EQ.  MAI  .OR.  I  .EQ.  MM2)  GO  TO  64 

62  EF11 1 = ( 3 . D0*RLI+2 . DO’RMUI ) * 1 . 2D1 *EX I *HAK/RHC 
GO  TO  65 

63  EF111=(3.D0*RL1+2.D0*RMU1 )*1 . 2D1*EX1*HAK/RHC 
GO  TO  65 

64  EF 1 1 1 M  3 . DO  *  RLC+ 2 . DO  *  RMUC  >  * 1 . 2D1 *  EXC  *  H AK/ RHC 

65  IF < I  .LT.  Ml  .OR.  I  .GT.  M2>  GO  TO  66 
IF< I  .GT.  Ml  .AND.  I  .LT.  M2)  GO  TO  67 
IF< I  .EQ.  Ml  .OR.  I  .EQ.  M2)  GO  TO  68 

66  EFI  J= ( 3 . D0*RLI+2 . D0*RMUI >  *1 . 2D1 *EX I ‘HAK/RHC 
GO  TO  70 

67  EFI  J= ( 3 . D0*RLl+2 . D0*RMU1 ) *1 . 2D1*EX1*HAK/RHC 
GO  TO  70 

68  EFI J=  <  3 . D0*RLC+2 . D0*RMUC ) *1 . 2D1 *EXC*HAK/RHC 

70  IF < I  .LT.  MM1  .OR.  I  .GT.  MA2)  GO  TO  71 
IF< I  .GT.  MM1  .AND.  I  .LT.  MA2 )  GO  TO  72 
IF( I  .EQ.  MM1  .OR.  I  .EQ.  MA2)  GO  TO  73 

71  EF211=(3.D0*RLI+2.D0*RMUI )*1 . 2D1*EXI*HAK/RHC 
GO  TO  74 

72  EF21 1 = ( 3 . D0*RLl+2 . D0*RMU1 ) * 1 . 2D1 *EX1 *HAK/RHC 
GO  TO  74 

73  EF211=< 3 . D0*RLC+2. DO'RMUC ) * 1 . 2D1 *EXC*HAK/RHC 

74  IF( I  .LT.  M112  .OR.  I  .GT.  M212)  GO  TO  75 
IF( I  .GT.  M112  .AND.  I  .LT.  M212)  GO  TO  76 
IF < I  .EQ.  M112  .OR.  I  .EQ.  M212)  GO  TO  77 

75  EF221=(3.D0*RLI+2.D0*RMUI )*3 . 2D1*EXI*HAK/RHC 
GO  TO  S4 

76  EF221  =  ( 3 . D0*RLl+2 . D0*RMU1 >  *1 . 2D1 *EX1 *  HAK/RHC 
GO  TO  S4 

77  EF221=( 3 . D0*RLC+2 . D0*RMUC ) * 1 . 2D1 *EXC*HAK/RHC 
GO  TO  S4 

53  EF121= <  3. D0*RL2+2 . D0*RMU2 ) • 1 . 2D1 *EX2*HAK/RHC 
EF111=( 3 . D0*RL2+2 . D0*RMU2 ) *1 . 2D1*EX2*HAK/RHC 
EFIJ= ( 3. D0*RL2+2 . D0*RMU2 ) *1 . 2D1*EX2*HAK/RHC 
EF211 = ( 3 . D0*RL2+2 . D0*RMU2 ) *1 . 2D1 *EX2*HAK/RHC 
EF221=( 3 . D0*RL2+2 . D0*RMD2 ) *  1 . 2D1*EX2*HAK/RHC 

54  IX=II+J-1 

IF< I  .LT.  IX)  GO  TO  30 
IF< I  .GT.  JJ)  GO  TO  31 

FX< I , J)=(EF211*T< I+l , J)-EF111*T< 1-1 , J) >/( 2.D0*XS< I , J) ) 

GO  TO  29 

30  FX< I , J )=(EF121*T< 1-2, J )-4. D0*EF111*T< 1-1 , J )+3.D0*EFIJ*T( I , J ) ) 
A/<2.D0*XS( I,J) ) 

GO  TO  29 

31  FX< I , J)=<-3.D0*EFIJ*T< I , J)+4 . D0*EF211*T< 1+1 , J )-EF221*T( 1+2, J) ) 
A/(2.D0*XS< I, J) ) 

29  CONTINUE 
28  CONTINUE 

DO  32  J=2,J3 
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DO  33  1=1,11 
IF( J  .LT.  NN1 >  GO  TO  83 
IF( J  .EQ.  NN1 >  GO  TO  84 
IF< J  .GE.  Nl>  GO  TO  8S 

83  ES121  =  <  3. D0*RLl+2 . D0*RMU1 ) *1 . 2D1*EX1*HAK/RHC 
GO  TO  86 

84  IF( I  .EQ.  Ml  .OR.  I  .EO.  M2)  GO  TO  87 
IF < I  .GT.  Ml  .AND.  I  .LT.  M2>  GO  TO  88 
ES121=( 3.D0*RLI+2.D0*RMUI ) *1 . 2D1*EXI*HAK/RHC 
GO  TO  86 

87  ES121=(3.D0*RLC+2.D0*RMUC)*1.2D1*EXC*HAK/RHC 
GO  TO  86 

88  ESI  21  =  <  3 .  D0*RLl+2 .  D0*RMlll )  *  1 . 2D1  *EX1  *HAK/RHC 
GO  TO  86 

85  ES121= ( 3 . D0*RL2+2 . D0*RMU2 ) *1 . 2D1 *EX2*HAK/RHC 
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THIS  SUBROUTINE  MAPS  TEMPERATURE  AND  ITS  GRADIENTS  IN 
THE  TEMPERATURE  FIELD  TO  THE  CORRESPONDING  POINTS  IN  THE 
STRESS  FIELD 


SUBROUTINE  MAP(TS, FXS,FYS,T,FX,FY, II , J1 ,LI1 ,KJ1 > 
IMPLICIT  REAL *8  (A-H.O-Z) 


T , FX , FY  =  TEMPERATURE  AND  ITS  GRADIENTS  IN  TEMPERATURE  FIELD 
TS.FXS.FYS  =  TEMPERATURE  AND  ITS  GRADIENTS  IN  STRESS  FIELD 


DIMENSION  T<I1,J1),FX<I1,J1),FY(I1,J1),TS<LI1,KJ1), 
AFXS(LI1 ,KJ1 ) , FYS(LI1 ,KJ1 > 


DO  1  J=1 ,KJ1 
DO  2  1=1, LI1 
TS< I , J )=0.D0 
FXS( I,J)=0.D0 
FYS( I , J)=0.D0 
CONTINUE 
CONTINUE 


DO  3  J=1,J1 
IT1=1 


DO  4  1=7,10 

TS< I , J) =T( IT1 , J ) 

FXS( I , J )=FX< IT1 , J ) 

FYS(I,J)=FY(IT1,J) 

IT1=IT1+1 


IT2=9 

DO  S  1=13, IS 
TS<  I , J)=T( IT2, J) 
FXS ( I , J ) =FX ( IT2 , J ) 
FYS< I , J ) =FY ( IT2 , J  > 
IT2=IT2+S 


IT3=22 

DO  6  1=16,19 
TS( I,J)=T< IT3,J) 
FXS< I , J  >=FX< IT3, J) 
FYS< I,J)=FY( IT3,J) 
IT3=IT3*3 


IT4=38 

DO  7  1=22,30 
TS< I , J )=T< IT4, J) 
FXS(I,J)=FX< IT4, J) 
FYS< I , J)=FY ( IT4, J ) 
IT4=IT4+3 


ITS=69 

DO  8  1=33,39 
TS< I , J ) =T< ITS, J  > 


-180- 


FXS< I , J)=FX( ITS, J) 
FYS< I ,  J ) =FY  < ITS , J ) 

8  ITS=ITS+3 
C 

ITS=91 

DO  9  1=40,42 
TS< I , J) =T< IT6, J  > 
FXS< 1 , J )=FX< IT6, J ) 
FYS( I , J) =FY ( IT6, J  > 

9  IT6=IT6+4 
C 

IT7=104 
DO  10  1=43,49 
TS< I , J)=T( IT7 , J  > 
FXS< I,J)=FX( IT7, J) 
FYS< I , J)=FY  < IT7, J> 

10  IT7=IT7+S 
C 

IT8=136 
DO  11  I=S0,S4 
TS(  I , J)=T( IT8, J ) 
FXS< I,J)=FX( IT8,J> 
FYS( I , J)=FY  < IT8, J) 

11  IT8=IT8+2 
C 

TS(11,J)=T<6,J) 
FXS(11,J)=FX(6,J) 
FYS ( 1 1 , J ) =FY  <  S , J ) 
TS(12,J)=T(8,J> 
FXS<12, J)=FX<8, J) 
FYS ( 12 , J ) =FY  <  8 , J ) 
TS ( 20 , J ) =T  <  33 , J ) 
FXS ( 20 , J ) =FX ( 33 , J ) 
FYS<  20 , J ) =FY  <  33 , J ) 
TS(21,J)=T<3S,J) 
FXS ( 21 , J  >  =FX ( 35 , J  > 
FYS  <  21 ,  J )  =FY  <  35 ,  J ) 
TS<31,J)=T<64,J> 
FXS ( 31 , J ) =FX ( 64 , J ) 
FYS ( 31 , J  >  =FY  <  64 , J  > 
TS(32, J)=T<66, J) 
FXS( 32, J)=FX(66, J) 
FYS( 32 , J ) =FY  <  66 , J ) 
3  CONTINUE 
C 

RETURN 

END 


MAIN  PROGRAM 

STRESS  FIELD  OF  A  LAYERED  MEDIUM  WITH  A  CAVITY,  THE  TOP 
EDGE  OF  THE  CAVITY  IS  AT  THE  INTERFACE 


C 
C 
C 
C 

IMPLICIT  REAL *8  <A-H,0-Z> 

C 

C  X  8.  Y  =  COORDINATES  IN  THE  PHYSICAL  PLANE 
C  A  =  MATRIX  TO  BE  SOLVED 

C  B  =  RIGHT  HAND  SIDE  OF  THE  ALGEBRAIC  EQUATIONS 

C  U  8.  V  =  DISPLACEMENTS  IN  X  AND  Y  DIRECTION,  RESPECTIVELY 

C  Sll ,S12,S22  =  STRESSES 

C  TS,FXS,FYS  =  TEMPERATURE  AND  ITS  GRADIENTS  FROM  TEMPERATURE  FIELD 
C 

DIMENSION  X< 67, 35 ), Y< 67, 35 >, A< 649740 ),B< 4420 >,TS< 67,35), 
AFXS(67,3S) ,FYS<67,3S) ,U<67,3S> ,V(67,3S) ,S11(67,3S) , 

AS12( 67, 35  > ,S22<  67 ,35) 

COMMON  /Zl/  RMU1 ,RMU2,RMUI ,RMUC,RMUIC,RMUI1 ,RMUI2,RMUC1 ,RMUC2 
COMMON  /Z2/  RL1,RL2,RLI ,RLC,RLIC,RLI1,RLI2,RLC1,RLC2 
COMMON  /Z3 /  DNi , DN2 , DN3 , DI , DJ , RHC 
C 

C  LI1  8.  KJ1  =  THE  TOTAL  NUMBER  OF  GRID  POINTS  IN  X  AND  Y  DIRECTION 

LI1=67 

KJ1=35 

Ll=LIl-2 

L2=LI1-1 

K1=KJ1-1 

C 

C  MB AND  =  HALF  BANDWIDTH 

C 

MBAND=K1* 2+6-1 
C 

C  NEQ  =  TOTAL  NUMBER  OF  EQUATIONS  TO  BE  SOLVED 
C 

NEQ=L1*K1*2 

C 

C  NTOT  =  TOTAL  DIMENSION  OF  "A"  VECTOR 
C 

NTOT=< 2*MBAND+1 >*NEQ 
JJ=K1*2 
C 

C  Ml  8.  M2  =  X  COORDINATES  OF  THE  CAVITY  CORNERS 
C 

Ml=19 

M2=31 

MR1=M1+1 

MR2=M2+1 

C 

C  N1  8>  N2  =  Y  COORDINATES  OF  THE  CAVITY  CORNERS 
C 


